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A STUDY AND EVALUATION OF IMAGE ANALYSIS TECHNIQUES 
APPLIED TO REMOTELY SENSED DATA 

By 

R. J. Atkinson, B,V. Dasarathy, M. Lybanon, H.K. Ramapriyan 


Abstract 


Several areas of concern to users of remotely sensed data, are covered 
briefly in this report. An analysis of phenomena causing nonlinearities in the 
transformation from Landsat multispectral scanner coordinates to ground coor- 
dinates is presented. Experimental results comparing rras errors at ground con- 
trol points indicate a slight improvement when a nonlinear (8-parameter) trans- 
formation is used instead of an affine (6-parameter) transformation. The 
improvements are expected to be more significant when the coordinates of the 
GCP’s are determined more accurately. Using a preliminary ground truth map 
of a test site in Alabama covering the Mobile Bay area and six Landsat images 
of the same scene, several classification methods are assessed. The similarity 
measures indicate a slight improvement over those with' a smaller test site in 
Tennessee used in an earlier report. However, due to the size of the present 
data set,registration has been a more difficult .problem, and with improved GCP 
selection, it is expected that the similarity measures will show further increases . 
A methodology has been developed for automatic change detection using classifi- 
cation/cluster maps. 

A sophisticated coding scheme is employed for generation of change depic- 
tion maps indicating specific types of changes such as from only the interior points 
of a given set of classes in map 1 to a set of classes in map 2 . Inter- and intra- 
seasonal data of the Mobile Bay test area have been compared to illustrate the 
method. A beginning has been made in the study of data compression by applying 
a Karhunen-Loeve transform technique to a small section of the test data set. 

The second part of this report provides a formal documentation of the 
several programs developed for the analysis and assessments presented here. 
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1. INTRODUCTION 


Earth related applications of space technology have been receiving consid- 
erable attention and support during the recent years. An important example in 
this area is the use of satellites such as Landsat and Slg^lab for remote sensing 
of the earth's resources. With the regular coverage of earth by Landsats 1 and 2, 
large quantities of data are collected and transmitted. A major concern for a 
user of such remotely sensed data is the accuracy of the results in terms of spatial 
registration and interpretation of multispectral signatures. A natural consequence 
of the repeated coverage of regions on earth by the Landsats is the ability to detect 
changes. An additional matter of interest to both the users and the agencies 
collecting and storing the data is the loss in accuracy that one might suffer if the 
data were compressed in order to achieve economies in transmission and archiving 
of the enormous volume of information. 

The Image Coding Panel of the National Aeronautics and Space Administration 
addresses the problem of assessing existing digital computer techniques covering 
the above aspects of image analysis. Reference [1] is an example of the assess- 
ment of classification techniques applicable to remotely sensed multispectral 
images . 

This report, a result of a three months' study* by Computer Sciences 
Corporation, is a continuation of the effort in [1] and covers a broader spectrum 
of analysis techniques. 


A large region, of approximately 1.4 million Landsat pixels, covering the 
Mobile Bay area in Alabama was used in the study as a test site. Six Landsat 
images of the same region were used from Computer Compatible Tapes (CCT's) 
corresponding to six different passes. This was a much larger data set than that 
considered in [1] which had only 52, 000 pixels. Therefore, the registration of 
the images to a common frame of reference was more difficult. The problems 
involved, mainly leading to a nonlinear (rather than linear) coordinate trans- 
formation, are discussed in Chapter 2. Also included in Chapter 2 is an approach 
to the implementation of nonlinear transformations on large data sets . Chapter 
3 gives a brief description of the various classification techniques that were applied 
to the data sets. 

The assessment methods employed here follow those in [1]. The classifi- 
cation maps are compared with a preliminary ground truth map (GTM) of the 
region. Due to the type of GTM available, certain special steps were necessary 


* Supported by NASA Contract NAS 8-32 107 . 
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in converting it into a digital format and preparing it for automatic point-by-point 
comparison with the classification maps. Chapter 4 presents the details of the 
preparation of the GTM and the results of the comparisons. 

A methodology for automatic change detection is developed in Chapter 5 
demonstrating the results of comparing inter- and intra- seasonal pairs of classi- 
fication maps . 

• A beginning has been made in the evaluation of data compression techniques . 
Application of a typical technique (Karhunen-Loeve Transform) to a small section 
of the Mobile Bay test site is illustrated in Chapter 6. 
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2. REGISTRATION AND GEOMETRIC CORRECTION 
2 . 1 Introduction 

In performing multitemporal classifications, change detection, and assess- 
ment of classification accuracies relative to the ground truth by a point-by-point 
comparison with a ground truth map (GTM), it is necessary to register the various 
scenes under consideration relative to a common basis. It is useful to employ 
the UTM (Universal Transverse Mercator) System as the basis relative to which 
all images are corrected. The GTM is generally not expected to have any dis- 
tortions other than a slight rotational error due to misorientation while digitizing. 
However, significant distortions are present in Landsat imagery and when several 
views of a given ground scene are to be handled it is possible that the distortions 
are different for different views. 


Experiments have shown that for small subsets of a Landsat frame (say, 
less than 500x500 pixels), an affine transformation of the form 


"I 

R 

= A 

n 

E 


Rq 

C 


N 


Co 

- ^ 




^ - 


where (R, C) and (E, N) are the pixel and UTM coordinates, respectively, provides 
an adequate model for the distortions and results in errors of less than one pixel. 
However, for larger subsets, such a transformation produces larger errors. 
Therefore, either several affine transformations with image segmentation or a 
single transformation with additional terms should be used. 

This chapter provides an analysis of the causes of the geometric distor- 
tion, showing why the linear model is not adequate. Also, experimental results 
are shown to illustrate the improvements in rms errors at a given set of Ground 
Control Points (GCP) when EN terms are included in the transformation. It is 
found, however, that the gain in accuracy is not sufficiently significant to justify 
application of the more complicated transformation in the present work. With 
greater accuracy in the determination of the GCP coordinates, the differences in 
rms errors between the linear and nonlinear models could be more pronounced. 

The results presented in the subsequent chapters of this report are based 
on affine transformations applied to all the data sets. But, a system of programs 
has been developed for appl5d.ng the nonlinear transformation so that when the 
GCP's are identified more accurately, the results can be updated with the new 
geometric corrections. The philosophy of these programs is slightly different 
from that of the affine transformation routines and is described briefly in Section 
2.5. 
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2.2 Effects Leading to Geometric Distortions of Im^es 
2.2.1 Effect of an Altitude Variation 

Some useful equations relating image coordinates to ground locations are 
given in Ref. 1. The coordinates (x,y) of a point in the image are given by 

X = (f tan -q) sin i|r 

y = (f tan 'q) sin (1) 

where y is along the (image of the) heading line, x is in the perpendicular direction, 
and the origin is at the image center (the image of the subpoint), (Eq. (1) actually ' 
defines the instantaneous field of view of the sensor. Effects such as the skew 
produced as a result of the finite time required for the Landsat multispectral 
scanner to form an image are not included.) The other quantities are defined 
as follows: f is. a scale factor, r| is the nadir angle subtended at the satellite by 
the great circle arc connecting the subpoint and the observed point and 1* is the 
azimuth of the observed point measured from the heading line (vertex at the 
subpoint). A spherical Earth is assumed, and the boresigjht direction of the scanner 
is assumed to be straight downward. 

The nadir angle "n is given by 


tan 


T1 = 


R sin 6 

R (1-cos 6) + H 


( 2 ) 


where R is the Earth's radius, H is the orbital altitude, and 6 is the geocentric 
angle subtended by the great circle arc connecting the subpoint and the observed 
point. Then the equation for position along a scan line, measured in the image, is 


X = astos ,3) 

E (1.00S 6) + H 

where R6 measures the distance on the ground from the subpoint to a point on the 
scan line. Here x=0 gives the center pixel of the line. 

The dependence of T| upon time is a property of the scanner mechanism. 
The question of interest here is: For a given ri (therefore a given x), how does 
(angular) position on the ground vary with changes in altitude H ? Differentiation 
of Eq. (3) produces 

^ „ sin 6 , . 

dH R (cos 6 - 1) + H cos 6 
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Letu=R6 denote a distance on the ground. Then the shift in ground position 
corresponding to a given pixel x, due to a variation in satellite altitude H, is 

^ ^ R sin (u/R) (5) 

dH R [cos (u/R) - 1]+H cos (u/R) 

Prom Eq. (5), for a 10 km change in altitude (a reasonable figure) the ground 
position corresponding to the end pixel on a line shifts by 1. 02 km, the equiva- 
lent of 17.9 pixels. However, u/R =0.0145 radians for this case (and smaller 
for other points), so small-angle approximations for the trigonometric functions 
are valid. If terms no higher than the first degree in u/R are retained, Eq. (5) 
reduces to 

du _ u (6) 

dH H 

Therefore, although this effect may be large it is approximately linear. So the 
linear transformation of Ref. (1) should account for it. 

For the dimension normal to scan lines, the analysis is different. The 
Landsat scanner sweeps out scan lines in groups of six simultaneously. Within 
each group of six scan lines the aspect angle effect prevails. However the ground 
position of each successive group of scan lines depends on the along-track velocity 
and the constant time interval . So for points more than a few scan lines apart in 
the along-track direction the latter effect predominates. 

From celestial mechanics it is known that (neglecting non-central gravita- 
tional terms and comparable effects) the orbital period is proportional to the 3/2 
power of the semimajor axis. For circular orbits, the velocity is constant and 
inversely proportional to the period. (Assuming a spherical earth, this is also 
true of the ground track velocity, which is of interest here.) So the average 
along-track distance v on the ground between scan lines is 

V = C (7) 

Since the time interval between scan lines (1/6 the time between scans of success- 
ive groups of six lines) is constant. In Eq. (7) C is a constant (a combination of 
several factors) and a is the orbit's semimajor axis - or radius, since the orbit 
is assumed circular. Therefore the fractional change in v due to a change in 
altitude is 


-3/2 da ^ -3/2 dH 

V a a 


( 8 ) 
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For a scan line increment nominally encompassing the same interval along track 
as the previous cross - track example, the shift due to a 10-km altitude variation 
is 190 m, the equivalent of 2.41 scan lines. Therefore an altitude change produces 
a smaller effect along-track than cross-track. Also, again the effect is a linear 
scale change. 

2.2.2 Effect Due to V arying Position Along Scan Line 

The nominal Landsat orbit has a ground track that repeats identically every 
eighteen days. However, because of small orbit variations the actual ground track 
may shift laterally, for several passes over the same region. Therefore a land- 
mark that is directly beneath the satellite's path during one pass may be off to the 
side during another . Because of this situation, it is of interest to consider the 
variation in the mapping from the ground to the image plane as a function of 
position along a scan line . 

To consider this effect, we return to Eq, (3). The quantity of interest is 
dx/du = (1/R) dx/d6 . From Eq. (3), it is given by 


^ j H-2(R+H) sin^ (u/2R) j 

du ^ I [2R sin2 {u/2R) + H]2 j 


For points within Landsat images the trigonometric terms are much smaller 
than the others. So, to a good approximation. 


dx _ f 
du ~ H 


(10a) 


X = (f/H) u = (constant) u (10b) 

In this approximation, the mapping from u to x is linear; projective effects do not 
cause any distortion. The error in this approximation is less than 1 percent; small, 
but not negligible. 

The above analysis should also apply, at least approximately, to a sensor 
boresi^t (pointing) error due to a nonzero satellite roll attitude. Therefore, the 
nonlinearity may be increased by some factor greater than unity. 

2.3 Discussion 

The distortion-producing effects on Landsat imagery of some physical 
processes that can be analyzed easily have been studied. It was found that the 
effects of orbital altitude variations and the projection into the image plane from 
various points along a scan line can be accounted for very nearly by a linear 
transformation, as was concluded in Ref. 1. 



However, it can be seen, that the approximations are not perfect. First the 
effects of altitude variation will be considered. It should be noted that Eq. (6), 
giving du/dH, is a small-angle approximation to the more accurate expression, 

Eq. (5). For the numerical example considered, the discrepancy between the 
two formulas amounts to slightly less than Im - about 1-1/3 percent of a pixel. 
Clearly the small -an^e approximation is valid in this case. Another-more 
serious -approximation has been made, however. The right-hand side of Eq. (6) 

(as well as Eq. (5)) depends on H. So the Taylor series giving the change in u 
induced by a change in H contains nonzero terms beyond first order. The amount 
of the nonlinearity is estimated by the next term, which has a magnitude of roughly 
0.2 pixel for the case used as an example. This is small, but possibly of border- 
line significance if one is attempting to register images to subpixel accuracy. 

When the same analysis is applied to dv/dH, ^venby Eq. (8), the magnitude 
of the nonlinearity is smaller. For the same numerical example considered pre- 
viously, the nonlinear correction amounts to a fraction of a meter. So it may 
reasonably be treated as negligible. 

The above analysis has considered the distortions in components (u,v) of 
locations on the Earth's surface mapped into given positions in the image, where 
u is in the scanning direction and v is parallel to the motion of the subsatellite 
point. Standard geographic coordinates used, such as the E and N of the TJTM 
system, are not parallel to u and v. So both E and N contain contributions from 
u, which the above analysis shows may lead to small but possibly sigtiificant 
nonlinearities in the mapping from the ground into the image. 

The other effect considered, the (projective) effect due to varying position 
along a scan line, does not involve changes in altitude. So Eq. (10a) actually 
describes a linear effect. But this equation is a small-angle approximation to the 
more accurate Eq. (9). Differences between the two equations give the sought- 
for nonlinearity. As has been pointed out, Eq. (10a) is a very good approximation, 
in error by only a fraction of a percent at the ends of a scan line . However, this 
translates into an .equivalent ground distance of about 200 m at these points, a 
significant nonlinearity equivalent to several pixels . And, as discussed above, 
the effect of this nonlinearity is felt in both E and N. 


It has been seen that, while the linear coordinate transformation is a very 
good approximation, the effects considered lead fco some possibly significant non- 
linearities. And only a few effects have been analyzed. It is possible that other 
causes may lead to even greater nonlinearities. Therefore it appears that a 
nonlinear transformation might be used with profit. 
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The simplest nonlinear transformation is obtained by adding'terms propor- 
tional to the product (EN) to both the x and y equations of the linear transformation 
derived in.Ref. 1. This simplicity is consistent with the relatively small magni- 
tude of the expected nonlinearities. Also, Eifman et al reported success in using 
such a transformation with Landsat data. [2] 


2.4 Examples 

Least-squares fits were made to several sets of data resulting from Landsat 
passes' over the Mobile Bay, Alabama area, using both a linear coordinate trans- 
formation and the nonlinear transformation described above . The linear trans- 
formation involved six parameters: the four elements of a matrix accounting for 
rotation, skew and scale factors, and the two components of an origin shift vector. 
The nonlinear transformation added two more parameters: coefficients of cross- 
product (E N) terms . The program used made several fits for each data set, 
successively discarding fit points possessing fit errors too large to pass a test. - 

The same 36 ground control points (GCPs) were used for each data set. 

Table 1 shows a sample of the results, for three data sets. RMS errors are- ' . 

given for fits involving all 36 GCPs and fits using 18 points (in some cases, linear' 
interpolation was used to obtain the latter) , Shown are errors computed at fit 
points only, and for all the points. 

The results show that the nonlinear transformation always produced smaller 
errors at fit points. (They can never be larger, since the nonlinear transformation 
was produced by adding terms to the linear transformation.) The situation is less 
clear-cut when errors computed for all 36 GCPs are compared for the 18-point 
fits. In order to clarify the situation, several steps may be useful: 

(1) Recheck the coordinates of all the GCPs. 

(2) Obtain better estimates of GCP image coordinates by first 
making enlargements (using digital resampling techniques) 
of regions surrounding the points, 

(3) Experiment with more data sets . 

(4) Experiment with other nonlinear transformations. 

Because all of these will involve significant effort and time, an estimate of the 
expected improvement over present results should be made and balanced against 
the expense required to achieve that improvement. 
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Table 2-l« ElVK Errors in Geographic Referencing Fits 
to Landsat Observations of Mobile Bay Area 


Date of 
Observation 

Number 
of GCPs* 

6-Parameter 
Linear Transformation 

8-Parameter 

Nonlinear Transformation 

Error at 

Pit Points All Points 

Error at 

Fit Points All Points 

10/17/72 

36 

1.7912 

1.7912 

1.7831 

1.7831 


18 

0.80430 

1.9095 

0.69475 

2.0202 

11/17/73 

36 

1.3308 

1.3308 

1.2320 

1.2320 


18 

0.73284 

1.4361 

0.57837 

1.4444 

12/5/73 

36 

1.3264 

1.3264 

1.2970 

1.2970 


18 

0.58774 

1.4614 

0.53989 

1.4512 


Rms errors are given in pixel units. 


*GCP: ground control point. 
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2.5 Geometric Correction 


In this section we shall discuss a method of appljdng a nonlinear geometric 
transformation to a large image. 

If the main memory of the computer is large enough to contain all the lines 
of the input image required to compute any one output line, then the implementation 
of the geometric correction is a simple matter. But, this is not the case in general. 
A method of handling large images using segmentation of the input image has been 
described in detail in [3] for the case of the affine transformation. This tech- 
nique is applicable to general nonlinear transformations, but the computations of 
the number of segments required and the portions of each output record that can 
be computed with the data from an input segment are somewhat cumbersome. 
Therefore, a slightly different approach is employed which will be described in 
the sequel. 

2.5.1 Statement of the Problem 


We are given a transformation of the form 


X 

y 



uv 


( 11 ) 


where X, Y are the coordinates of a. point in the input image reference system 
and U, V are those of the same point in the output coordinate system. The input 
image function (say classification numbers or reflected intensity values in a given 
spectral band) is defined for all integral values of X, Y satisfying 

1 ^X^X; 1sY:£Y (12) 

The output image function should be stored as (U - U + 1) records (lines) with 
(7- V + 1) samples in each record, the values being computed for all integral 
values of U, V satisfying 

U^U^U;V£V^V' (13) 

where U , U , V , V are the integral bounds on U and V such that (11) and (12) are 
satisfied. 

The following steps are involved in computing the output image function. 
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1. Determination, of U , U, V, V. 

2. Deciding whether it is convenient to implement the given trans- 
formation or a slight modification thereof. 

3. Computing the input image coordinates for each output point in 
a record. 

4. Finding the nearest integer values to the input coordinates and 
obtaining the output image function by a suitable interpolation method. 

2.5.2 Determination of Bounds on the Output Image 

The values of U , tJ , V , V should be found such that, for all (U, V) 
satifying (13), the solution (X, Y) to equation (11) satisfies (12). 

If 0 ! = P = 0 the values of U , U, V , V are simply found by using 


u 

1! 

> 

X -x„l 

V 


1 

o 

r 

1 


(14) 


and finding (U, V) corresponding to the four vertices of the rectangle defined by 

( 12 ). 


However, equation (11) results in quadratic equations for U and V in terms 
of X and Y. The equation for U is: 


(aj^lp-agi Q!)u 2+ [aj^322“aj^2 ^21'*' Q;(Y-Yq)- P(X-Xq)']u 

+ [ai2(Y-Yo)-a22(X-Xo)]=0 (15) 


This 3 delds two solutions for U. 

A choice is made of the sign to be used in the formula for the solution of 
the quadratic equation (15) based on the absolute value of the solutions for the case 
(X, Y) = (l, 1). The sign yielding the smaller absolute value is chosen for con- 
venience and is used henceforth. Given a U, the solution for V is unique.. 

Due to the nonlinear nature of the transformation, the four vertices of the 
rectangle defined by (12) do not necessarily yield the extremal values of U and V . 
However, an examination of the quadratic curves represented by the equation (11) 
in the (U, V) domain for various values of (X, Y) reveals that the sides of the 
rectangle defined by (12) should map into curves bounding the output image. 
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Therefore, U, U, V, V are obtained by finding the extremal values of U and V 
over the sets 


x = l ; 

1 ^ Y £ Y 

X =X; 

1 ^ y £ Y 

1 ^X s 

X ; Y = 1 

1 ^X ^ 

X; Y = Y 


2.5.3 Modification of the Transformation 


( 16 ) 


It is sometimes convenient to implement a modified transformation and 
generate an image which is equivalent to the desired output image from the point 
of view of resampling (that is, integer coordinates in the generated output image 
correspond to integer coordinates in the image desired). The permissible modifi- 
cations are interchanges of the IT and V axes and/or changes in the signs of U and/ ‘ 
or V. 


The maximum number of input records required to compute an output record 
is given by 

%=Max[^ja;^2‘^‘^^l| I ai2 + «U | (V-V)J 

if the given transformation is implemented. If U and V axes are interchanged, 
this changes to 

Rg = Max y ^ I (U -U) I &ii + 0! V j (U -U)J 

If R£ >_fhjen U and V need not be interchanged. Otherwise, the columns 
of A, (U, V) and (U, V) are interchanged. 

Also, it is desirable to have the output record and sample numbers increase 
with respect to those in the input image. Since a and ^ are generally much smaller 
than the elements of A, a^i S 0 and S 22 ^ 0 imply that the above requirement is 
met. If an < 0, then the first column of A, a, P, U, U are negated and U, U 
are interchanged. If 022 < second column of A, a , p, V, V are 

negated and V, T are interchanged. 

2.5.4 Computation of Input Image Coordinates 

For each line in the output image (that is, for a given U) there are V- V + 1 
output pixel values to be computed. The (X, Y) coordinates for each of these 
points are computed and stored in two arrays. Also, some of the computation is 
saved by checking whether the computed value of X is between 1 and X. If not, 
the Y value is set to 0 without any computation. 
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2.5.5 Computation of Output Pixel Values 

Suppose |Xj I i = l, 2, N} and{Yj|i = l, 2, n} are the coordinate 

arrays with (Xj[, Yj) cooresponding to (U, i=V-l) andN=V-V+l. Also, when- 
ever 1 >Xj or Xj > X, Yj=0. Then, the input records needed to compute the 
current record of output range from throu^ R 2 where 

Rl = Min jXjj 1 < i s N ; Yj 0^ 

R2 = Max|Xj|l ^i ^N;Yj?^of. 

Now, if the storage available for the input image can contain R 2 - R^ + 1 
records, it is possible to compute the entire output record by reading the appro- 
priate input records into storage. However, if only apart of the records needed 
can be held in core at a time, the data handling becomes more complicated. The 
method employed here is iterative. 

Assume that records through S 2 are available in core at a given stage 
and Rl through R 2 are needed for computing the remaining part of the current . 
record.’ Then, all the output pixels that can be computed using Sj through S 2 are 
computed, a code array indicating the computed pixels is updated, as many as 
possible of the new input records in the range R^ through R 2 are read into core 
and Si, Sg are modified to reflect the present state of the storage. This procedure 
is repeated until all the output pixels of the current record have b een c omputed. 

The above method works employing a circularly addressed input buffer to 
avoid unnecessary movement of data within the core storage and assumes that the 
input image is available on a direct access device so that the updating of storage 
can be performed using both forward and backward reading of input records . 

2.5.6 Experimental Comparison of Transformation 

An experiment was performed using both the affine and the above 8-parameter 
transformation on a six class classification map of the Mobile Bay, Alabama region 
obtained by using a Linear Sequential classifier. The transformations used here 
were both obtained such that the mean squared errors over the same set of ground 
control points were minimized. The resulting output images were overlaid and 
their joint histogram was found. This is shown in Figure 2.1. The "sindlarity 
measure” between the two images, defined as in M is seen to be 84.6 percent. 

Also, as ejqjected, most of the deviations occur at boundaries between classes. 

In fact the differences occuring at locations Interior to the classes in the output 
of the affine transformation amount only to 2.79 percent. 
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3. 


PKOGRAM DESCRIPTIONS 


3 . 1 Introduction 

In this study, three classification techniques in addition to those described ' 
previously [l] were employed. These additional techniques are supervised, in 
that a set of training samples, whose classification is known, is input to the sys- 
tem. These samples are then used to evaluate the required parameters of the 
classification algorithm, such as the Gaussian parameters of the distributions or 
the coefficients of linear discriminant functions, for use in a table look-up classi- 
fication scheme. 

The next step is the table formulation phase, which consists of creating 
the table or array of class labels corresponding to the set of all the independent 
feature vectors expected in the data environment. In the approach considered 
here, this effort is limited to generating the cluster or class labels corresponding 
to a smaller set of prototypes only instead of all possible feature measurements. 
The prototypes are defined here as the centroids of the contents of all the occupied 
cells resulting from a multidimensional histogram analysis. When being deployed 
with the HINDU [2] system, thei’e is no additional effort needed for identifying these 
prototypes as these are already available being the output of the histogram gener- 
ator. 


The classification, or table lookup, phase requires the use of the incoming 
measurement vector to locate the proper element of the table, which contains the 
class number to which the pixel is assigned. 

In the approach developed here for use with the HINDU system, each input 
sample has already been flagged during the histogram analysis with an index- 
number which identifies its address on an address array, this address denoting 
the cell corresponding to the input sample. This information being available on 
external storage is read in (instead of the feature vectors themselves) during 
the table lookup phase and the entries in this array are replaced by the corres- 
ponding entry of the label array or table. Hence, no additional searching of the 
array is involved in this phase (equivalent effort having been expended earlier in 
the learning phase) . 

3 . 2 Resource Requirements 

The computer resources employed by these classifiers are similar to those 
given previously for the HINDU system, viz . a core memory requirement of 
15 OK bytes, two tape drives for input and output, and external storage for the 
histogram cell addresses. An additional input requirement is a set of labeled 
training samples . 
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3.3 Supervised Nearest-Neighbor Establishing Histogram Approach (SNEHA) [3] 
3.3.1 Analysis Process 

In this method, each class is represented not by a single description (e. g, 
its centroid), but by a set of descriptors. Here, the centroids of all the histogram 
cells within each class are identified and designated as the descriptors of the 
corresponding clusters. The prototypes are then classified on the basis of the 
label of the nearest descriptor (cell centroid). 

The steps in the analysis are as follows: 

(i) Training Set Definition: This phase is common to all supervised 
classification techniques and accordingly details of these efforts 
are omitted here, except to state that training samples corres- 
ponding to all of the pattern classes expected in the data environ- 
ment are selected through ground truth and/or photo -interpretation 
techniques . These samples along with their labels form one of 
the major inputs to the ensuing phases. 

(ii) Prototypes Identification: Through a multidimensional histogram . 
analysis of the entire data set, the centroids of all the occupied 
histogram cells are extracted. These centroids then represent 

a set of prototypes of the given unclassified data set. These 
prototypes are treated as a pseudo sample set to be classified 
according to the nearest neighbor rule. Further, each input 
sample is also identified during this process by the relative 
address of the multidimensional histogram cell occupied by the 
sample. This information is stored for all samples externally 
for later retrieval. Here, the computational effort is a function 
of several variables: The dimensionality and size of the data 
set, the grid size of the histogram. While the first two are 
dictated by the data set, the latter is an option available to' the 
user. This has been discussed at length in the HINDU System 
wherein a similar histogram analysis is performed prior to 
unsupervised learning. The grid size is therefore chosen in 
accordance with the guidelines and empirical relationships de- 
veloped for the HINDU system. It is to be noted that in general, 
the computational effort increases inversely with grid size and, 
hence, a trade off study between computational effort and accuracy 
may be necessary. 

(iii) Label Table Formulation: The training sample set along with its 
labels forms the data base for the nearest neighbor classifier. 

The prototypes, i.e. , the centroids of all the occupied histogram 
cells are classified on the basis of nearness to one or the other of 
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these training samples, and the corresponding labels are stored 
in the form of a table or array of labels to be used in the table 
lookup phase. The computation involved in this phase is pro- 
portional to the size of the training sample set as well as the 
size of the prototype set, but is independent of the size of the 
total data set to be labeled. This independence is the key to the 
success of this approach in the processing of large data sets . 

(iv) Table Lookup and Labeling; In this phase, the relative address 
of each input sample is retrieved from the external memory 
(wherein it was stored during the histogram analysis process) 
and the corresponding entry in the label table is looked up to 
identify the label of the input sample. Repetition of this opera- 
tion for the entire data set leads to the output label set corres- 
ponding to the input unclassified sample set, 

3.4 Parametric Recognitionimparting Trained Identification (PRITI) [4] 

3.4.1 Analysis Process 

The classical maximum likelihood method is used to classify the centroids 
of the occupied histogram cells. The steps in the analysis are as follows: 

(i) Training Set Definition: This phase, being common to all supervised 
classification schemes, is not discussed here at length. The approach 
is, of course, based on the availability of well defined ground truth 
(i.e. , labeled training sample set) and a knowledge of the probabi- 
listic description of the pattern classes (i.e. , distributions under- 
lying the data). These form the input to the next phase. 

(ii) Parameters Estimation: This is the so-called learning phase, 
wherein the Gaussian parameters of the distributions underlying 
the different pattern classes are estimated using the labeled 
training sample data. 

(iii) Prototype Identification : The input data set consisting of all the 
unclassified samples is processed through a multidimensional 
histogram analysis package to identify the centroids of all the 
occupied histogram cells. These centroids are then considered 
as prototypes of the incoming data set to be classified by the 
classical parametric recognition method employing the maximum 
likelihood classifier. During this histogram generation, the 
relative address of the histogram cell corresponding to each in- 
put sample is developed and stored externally for future retrieval. 

The computations involved in this phase are proportional to the 
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size of the data set and represent the major part of the total 
computational expense. However, the histogram grid size has 
an equally significant effect on this computational load. This 
size is presently dictated by certain semi-empirical considera- 
tions developed earlier in the HINDU system of unsupervised 
learning. 

(iv) Table Formulation; The prototypes of the data set as discerned 
I by histogram analysis are classified here using the classical 

maximum likelihood classifier with the estimated parameter 
values. The labels corresponding to the prototypes are stored 
as a table or array of labels for the ensuing table lookup phase. 

The computations involved in this phase are proportional to the 
number of prototypes (and hence the grid size and spread in the 
sample set), square of the dimensionality of the data, and the 
number of pattern classes. 

(v) Table Lookup and Labeling! Here the relative address of each 
sample in the histogram space is recalled from external memory 
(where it was stored earlier in the histogram generation phase) 
and used to lookup the table entry to identify the label of the in- 
put sample. This process of table lookup is repeated for all the 
input unclassified samples to derive the output label set. 

3.5 Piecewise Linear Classifier 
3.5.1 Analysis Process 

The occupied cells are classified on the basis of their centroid positions 
by a set of discriminant hyperplanes. The parameters of the discriminants are 
determined by the set of training samples . Otherwise, the procedure is identical 
to that described in the preceding section regarding parametric recognition. The 
positions of the discriminant hyperplanes are computed by the Discriminant 
Hyperplane Abstracting Residuals Minimization Algorithm (DHARMA).[5J which has 
been described in conjunction with the HINDU system. Consequently, the analy- 
sis process is not repeated in detail here. 

3 . 6 Performance Assessments 

The CPU time required is divided between the histogram analysis time and 
the table creation and lookup time. Typically, processing a four-dimensional 
data set consisting of a quarter million pixels required 84 seconds for identifying 
the significant clusters in terms of their centroids . 

The CPU time required for creating the table containing the class numbers 
of the prototype cells is on the order of a few seconds, with a maximum of 
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approximately three seconds for the piecewise linear method. This is due to the 
longer time required to determine the linear discriminant functions, when com- 
pared to' the other methods. 

The table lookup phase is very rapid. The computational demands of this 
phase, in terms of CPU time, are only nominal with most of the effort being 
I/O operation (which is proportional to the number of records or scan lines of 
data being processed). Even compared to other table lookup approaches, the 
savings achieved here is significant as the number of prototypes which are 
classified by the classifier is still much smaller than the number of nonidentical 
feature vectors one comes across in a large Landsat scene. Approximately 1 1/2 
seconds of CPU time is expended in looking up the table to identify the class 
labels of the quarter million pixels . 
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4. CLASSIFICATION TECHNIQUE ASSESSMENT 


4-1 Introduction 

The approach used in this report for assessing the classification methods 
closely parallels that of [1]. The data sets used here are much larger in size, 
with approximately 1.4 x 10® pixels compared to 51000 in [1], The details of 
the data sets are given in the next section . 

The ground truth map (GTM), supplied by the U. S. Geological Survey* is 
a second level classification map in which each region is denoted by a two digit 
number. For the purposes of the present work this map was converted into digital 
format with only the first level classifications . The details of rendering the GTM 
to a digital form are presented in Section 4.3. 

Various tests have been applied to the classification maps described in 
Chapter 3 . The similarity measures relative to the GTM, typical contingency 
matrices, tests and inventory comparisons are shown in Section 4.4. 

The description of the data sets given below include tape numbers containing 
the data. These refer to the IBM 360 tape library in Building 4708, Marshall 
Space Flight Center, Huntsville, Alabama, and are provided only for the sake of 
completeness of the record. 

4.2 Description of Data Sets 

The data sets used here were obtained from the computer compatible tapes 
of the following six Landsat Multispectral Scanner (MSS) images , each with four 
bands of information: 

Table 4.1. Landsat Scenes of Interest 



ID 

Date 


1086-15562 

Oct. 17, 1972 


1482-15541 

Nov. 17, 1973 


1500-15535 

Dec. 5, 1973 

A0691 

1626-15510 

April 10, 1974 

A0484 

1698-15490 

June 21, 1974 

A1310 

1896-15420 

Jan. 5, 1975 


♦This map was preliminary and had a disclaimer indicating that it was not 
ready for public release. 
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The geographic region covered by the data sets is the Mobile Bay area in 
Southern Alabama. 

From each of the Landsat images, a 1200 by 1200 image was extracted to 
cover this region. Also, the "synthetic" pixels (SE) were removed after the ex- 
traction resulting in slight reduction in the number of pixels per line of the images. 
The coordinates of the top left corners of the data extracted relative to the 
corresponding Landsat images (before removal of SP's) and the image sizes (after 
removal of SP's) are shown in Table 4.2, 

These data sets were classified using various classification methods . Both 
the data sets and the classification maps were geometrically corrected to the 
UTM coordinate system using the affine transformations defined by the parameters 
shown in Table 4.3, and using a scale factor of 20 pixels per kilometer. Also, 
to register the geometrically corrected images with respect to the GTM prepared 
as described in Section 4.3, it was necessary to extract subsets of the images 
starting at the locations shown in Table 4.4. Finally, the registered images 
were all 1624 lines by 866 pixels in size. 

The tape numbers for the geometrically corrected data and the classification 
maps (GC) and the corresponding 1624 x 866 data sets overlaying the GTM (OV) 
are shown in Table 4.5. 

4.3 Preparation of the GTM 

The goal in preparing the GTM is to generate a digital data set consisting 
of the class numbers at each point of a rectangular array matching the UTM 
coordinates with 20 pixels /km (That is, each sample representing a 50 x 50 m^ 
area). Table 4.6 shows the various stages involved in the preparation of the ‘ 

GTM. In the following description GTMi will be used to denote the ith stage 
GTM. 


A preliminary map, GTMO, was supplied by J. R. Anderson, U. S. 
Geological Survey (with a disclaimer indicating that it was not ready for public 
release) . It consisted of interclass boundaries with the class identification 
numbers (GIN) written inside the respective regions. Being a level II classifi- 
cation map it was very detailed. Each GIN was a two digit number representing 
the level I land use class and the level II subclass . 

The portion of interest in GTMO was photographically copied and the 
annotations were masked out manually. The resulting image was photographically 
reduced to 4" x 3" (GTMl) for digitization on a microdensitometer. The image 
was digitized at a scanning interval of 25 microns to ensure that the boundaries 
were recorded properly. The part of interest in the map, designated as GTM2, 
has 4000 X 2100 pixels. 
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Table 4.2. Subsets of Landsat Scenes 


Date 

Initial 

Row 

No. of 
Rows . 

Initial 

Column 

No. of 
Columns 

Feature Vector 
Tape Numbers 

Bytes 

Floating Point 

Oct. 17’, 72 

156 

1200 

551 

1194 

A0923 

A0198 

Nov. 17, 73 

7 

1200 

665 

1193 

A1048 

A0787 

Dec. 5, 73 

62 

1200 

751 

1194 

A0052 

A0454 

Apr. 10, 74 

84 

1200 

788 

1194 

A1119 

A1235 

Jun. 21, 74 

148 

1200 

900 

1193 

A0872 

A0705 

Jan. 5, 75 

7 

1200 

836 

1193 

A0906 

A1347 


Table 4.3. Affine Transformation Parameters 


Date 

No. of 
GCP’s 

PARAMETERS 

RMS 

Error 

I— 1 

^12 

^1 

^22 

Rq 

Co 

Oct. 17, 72 

18 

-2.4037 

-12.2023 

16.6449 

-4.2916 

42758.3 

8863.9 

0.-80 

Nov. 17,73 

20 

-2.3846 

-12.2867 

16.6352 

-4.2836 

43098.4 

8838.8 

0,74 

Dec. 5, 73 

25 

-2.3416 

-12.3164 

16.6034 

-4.2711 

43192.2 

8798.4 

0.88 

Apr. 10, 74 

17 

-2.3816 

-12 1 1802 

16.6193 

-4. 1945 

42790.4 

8541.5 

0.83 

Jun. 21, 74 

24 

-2.4133 

-12^2759 

16.8070 

-4.4704 

43138.9 

9376.3 

0.79 

Jan. 5, 75 

24 

-2.3975 

-12.3053 

16.5819 

-4.2947 

43199.3 

8939.1 

0.78 





















Table 4.4. Information for Registration w. r. t. GTM 


Date 

Coordinates Matching (1, 1) on GTM* 

Row 

Column 

Oct. 17, 72 

248 

605 

Nov. 17, 73 

336 

577 

Dec. 5, 73 

345 

564 

April 10, 74 

430 

553 

June 21, 74 

431 

514 

Jan. 5, 75 

396 

615 


*The GTM here refers to GTMIO described in Section 
4.3 (Tape Number A0944) 




Table 4.5, Registered Data Sets 


Date 

Tape Numbers 

Raw Data 

HINDUCM 

SNEHACM 



PRITICM 

LCM 

PCM 

ETCM 

GC 

GC 

OV 

GC 

OV 

GC 

OV 

GC 

OV 

GC 

OV 

GC 

OV 

Oct, 17, 72 

AO 620 

A0869 

A0969 

A0557 

A0912 

A0575 

A0147 

A0528 

A0517 

A0390 

A0413 


- 

Nov. 17, 73 

A0671 

A0784 

A1209 

A0335 

A0492 

A0894 

A1437 

A0219 

A0202 

A0085 

A1498 

- 

- 

Dec. 5, 73 

A0063 

A0286 

A0200 

A0676 

A0633 

A0820 

A0822 

A0798 

A0790 

A0865 

AO 831 

A1152 

A1153 

Apr. 10, 74 

A1325 

A0966 

A0197 

A0259 

A0495 

A0209 

A0438 

A1434 

A1423 

A1222 

A1352 

- 

- 

June 21, 74 

A1287 

A0868 

A0846 

A0626- 

AOS 74 

A0440 

A0331 

A1457 

A1464 

A0573 

A0340 


- 

Jan. 5, 75 

AO 828 

A0109 

A0016 

A0279 

A0287 

A1388 

A1368 

A1315 

A0535 

A0498 

A0402 


- 







Now,- GTM2 had boundary lines with thickness greater than one pixei. (While 
thinner boundary lines could have been obtained by digitizing at SOM" resolution, 
undesirable discontinuities would have resulted in that case) . The boundary lines 
in GTM2 were thinned by a "peeling" algorithm which removes outer layers of 
"thick" lines while ensuring that connectivities are preserved. The resulting 
image, GTM3, was stored using the "scan line intersection code (SLIC)" wherein 
only the column coordinates of the intersections of each scan line with the boundary 
image are recorded. 

The information in GTM3 was converted into a Region Identification Map 
GTM4 wherein unique numbers (RIN) were used to identify each connected region, 
the boundaries being denoted by 0. A table of correspondences was established 
manually between the CIN’s of GTMO and RIN's of GTM4. It was found, however, 
that there were some small discontinuities in the digitized boundaries which 
caused some large distinct regions to be merged. Therefore, the discontinuities 
in GTM3 were automatically sensed and patched by generating straight line seg- 
ments. The result, GTM5, was converted into a Region Identification Map GTM6. 
The regions which were in error in GTM4 were now corrected, as far as possible 
by a new table of BIN to CIN correspondences. The two correspondence tables 
were used on GTM4 and GTM6 to get a map GTM7, storing the CIN's at all the 
points in the properly identified regions, 0*s at the boundary points and -1 at all 
the points in the regions still in doubt. The percentage of -I's in GTM7 is 2, 8. 
These points represent locations in the image which, due to errors in GTMO, 
indicated nonunique or unavailable assignments of CIN's. The corrections to the 
GTMO were received, but not in time to be incorporated into the present work. 

Next, GTM7 was corrected to UTM coordinates using an affine transfor- 
mation determined to minimize the mean squared error over a set of ground 
control points. The resulting images GTM8 was arranged to have 20 pixels /km 
in both coordinate directions. The GTM, at this stage, was still a level II 
classification map. For the present analysis, a level I map was required. There- 
fore, GTM8 was modified by table look-up to GTM9 with only 6 classes (1-Urban, 

2 -Agriculture, 3-Forest, 4-Water, 5-Wet land, 6-Vacant, 0-Boundary, Unknown 
or exterior of the image region after geometric correction). Next, each of O's 
in GTM9 was replaced by CIN occuring most frequently in the 3x3 region centered 
at it. Thus, O's representing the boundaries between subclasses of level I classes 
were eliminated while retaining the O's representing the unknown and exterior 
regions. The image, GTMIO, is shown in Figure 4.1. 

4.4 Comparisons With the GTM 

The map GTMIO was used as the standard for evaluating the classification/ 
cluster maps produced by the various methods. Only the points in GTMIO with 
class numbers 1 through 6 were used for determining the similarity measures . 

The following definitions of similarity measures parallel those in [l] . 
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Table 4.6. GTM in Various Stages of Preparation 


Stage 

Description 

0 

Preliminary map (with disclaimer) supplied by 
U. S. Geological Survey. 

1 

Annotations removed and photographically 
reduced . 

2 

Digitized at 25 M- scanning interval. 

3 

Boundary lines thinned. 

4 

Regions identified by unique numbers (RIN) 
based on connectivity. 

5 

Discontinuities in boundary lines patched. 

6 

Regions identified by unique numbers (RIN) 
based on connectivity. 

7 

RIN's replaced by CIN's using a manually 
determined correspondence table with regions 
in doubt flagged with CIN = -1. 

8 

Geometrically corrected to the orientation of 
UTM coordinates with 20 pixels Ani. 

9 

Reduced to Level I map by modifying the 
CIN's. 

10 

Boundary points replaced by CIN's occurring 
the largest number of times in their immediate 
neighborhoods . 
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Figure 4.1, Mobile Bay GTM 
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All the supervised classification maps (CM) have 6 classes - Urban, Agri- 
culture, Forest, Water, Wet land and Vacant - with the exception of ELLTABCM 
which has a seventh "reject" class. 

With A defined as the joint histogram between the GTM (rows) and a CM 
(columns), the similarity measure is defined as 

Tr (A) / S (A) 

6 

where Tr (A) = E a,-. 

i=l 

and S (A) = Sum of all elements of A, 

In the case of unsupervised CM's, the columns of A are interchanged and/or 
merged according to a reassignment rule and then the above definition is used to 
find the similarity measure. 

Due to uncertainties in registration and determination of geometric correction 
transformations, it is likely that the locations of the boundary pixels between 
classes in the GTM are unreliable. Therefore, the errors at the boundary and 
interior pixels are reckoned separately and similarity measures considering 
only the interior pixels are shown in addition to the "normal" similarity measures 
which consider the entire map. Also, "inventory" similarity measures are shown 
in each case indicating the closeness of the estimates of percentage occupancies 
of each Land Use Class in the CM's and the GTM, As in [l] , this similarity 
measure is defined as [1 ~ ElPii -P 2 il /2E Pj^j] 100 percent. In the case of the 
unsupervised maps the similarity measure is based on the same measurements 
as for "normal". 

In the following subsections, the CM's and the similarity matrices are pre- 
sented for each of the algorithms discussed in Chapter 3. With the exception of 
ELLTAB, the classification methods have been applied to all six data sets 
described in Section 4.2. However, in each case, only one typical CM and a 
similarity matrix are shown, with the similarity measures tabulated for all six 
cases. 

4.4.1 HINDU Classifier 

This histogram dependent unsupervised method generated classification 
maps with different numbers of clusters for the six data sets . The map for the 
January 1975 data set, consisting of 8 classes, is shown in Figure 4.2. The 
joint histogram of GTM v/s HINDUCM for January 1975 is shown in Table 4.7. 

This being an unsupervised map, the class numbers were assigned to numbers 
1 through 6 of the GTM classes and the similarity measures were evaluated. The 
similarity measures for the six HINDUCM's are shown in Table 4.8. 
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Figui'e 4.2. 


Mobile Bay HINDUCM (Jan. 1975) 
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Table 4.7. Joint Histogram Between GTM and HINDU Cluster Map 


class.no. 

1 . 

_ . 2 __ 

_3.58U 

. 3 

2679 

4 _ 

378 

5 

6 

7 

8 _ 


723A6 

281 A 

1567 

237 

5C93 

2 

81189 . 

_81745 

2G3 

31 __ 

6634 _ 

845 

1 140 1 

1742 

3 

46A658 

47892 

3on 

367 

2 31 5 

3562 

16 32 

16 52 

4 

4439 

16 31 19C25C 

J84327 

62 

12020 

^ ^ 

10 

L O y C 

1230 

_ 5. 

78792 

— -3436 

3656^ 

5C2 . 

42_ 

8899 

259 

297 

6 

4656 

_ -2CC6 

40C 

84 

.. 33 

10 .0 „ 

100 

^1912 


CT M V/S HlN OUCPfOl CS7SI 


Table 4.8. Similarity Measures of HINDUCM's w. r. t. GTM 


Date 

No. of Classes 

Similarity Measures (%) 

Normal 

Boundaries 

Ignored 

Inventory 



Oct. 17, 72 

9 

67.42 

69.90 

76.15 

Nov. 17, 73 

8 

70.96 

73.54 

90.45 

Dec. 5, 73 

9 

73.37 

76.08 

83.82 

Apr. 10, 74 

8 

69.67 

72.39 

78.80 

June 21, 74 

10 

69.06 

71.65 

83.63 

Jan. 5, 75 

8 

71.36 

74.06 

85.78 
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4.4.2 Table Look-up (ELLTAB) 

The classification map using this table look-up implementation of the 
maximum likelihood method for the December 1973 data set is shown in Figure 
4.3. It has 7 classes, the last being a "reject" class. The joint histogram of 
this map with the GTM is shown in Table 4.9. This method was used for the 
classification of only the December 1973 data. The similarity measures in this 
case were found to be: 

Normal: 73.48% 

Boundaries Ignored: 76.51% 

Inventory: 98.28% 

4.4.3 Linear Sequential 

The output, LCM, of this classifier is shown in Figure 4.4 for the April 
1974 data set. The joint histogram with GTM is shown in Table 4. 10, and the 
similarity measures are given in Table 4. 11 for all the six data sets. 

4.4.4 Piecewise Linear (Histogram based) 

This supervised classifier is implemented using a table look-up approach 
wherein, for a given grid size, the histogram cells corresponding to each of the 
feature vectors in the data set are identified first, a table of classifications of 
the occupied histogram cells is formed using training samples and a piecewise 
linear classification rule and, finally, the entire data set is classified by table 
look-up. The result, PCM, of applying this method to the October 1972 data 
set is shown in Figure 4.5. The corresponding joint histogram is shown in Table 
4.12 and the similarity measures are given in Table 4. 13. 

4.4.5 Nearest Neighbor (Histogram based) 

This is a supervised classifier similar in implementation to the piecewise 
linear classifier above. The table of classifications of the multidimensional 
histogram cells is produced using training samples and a nearest neighbor rule 
of class assignment. The result, SNEHACM, of this method is shown in Figure 
4. 6 for the November 1973 data set. The joint histogram and the similarity 
measures are shown in Tables 4. 14 and 4. 15 respectively. 

4.4.6 Maximum Likelihood (Histogram based) 

This classifier differs from the conventional maximum likelihood and 
ELLTAB in that it uses a multidimensional histogram with a specified grid size 
as in the case of PCM and SNEHACM. The resulting histogram cells are classified 
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Figure 4.3. Mobile Bay ELLTABCM (Dec. 1973) 
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Table 4-9. Joint Histogram Between GTM and ELLTAB Classification Map 
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Figure 4.4. Mobile Bay LCM (April 1974) 
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Table 4.10. Joint Histogram Between GTM and Linear 
Sequential Classification Map 


CLASS NO. 

1 

2 

3 

4 

5 

6 

1 

43249 

47676 

21100 

2129 

4462 

1411 

2 

23972 

133C34 

22448 

31 

3445 

960 

3 

41336 

1C5547 

359656 

564 

14984 

291 

4 

1183 

564 

1199 

386163 

5664 

1353 

5 

5991 

3239 

47f75 

33C9 

3551 5 

92 

6 

1957 

3526 

2589 

521 

622 

564 


GTM V/S LCMIOAICIA) 


Table 4.11. Similarity Measures of LCM's w, r. t. GTM 


Date 

Similarity Measures (%) 

Normal 

Boundaries 

Ignored 

Inventory 

Oct. 17, 72 

65.77 

68.57 

85.75 

Nov. 17, 73 

72.00 

74.95 

96.99 

Dec, 5, 73 

73.85 

75.81 

97.22 

Apr. 10, 74 

72.15 

75.14 

91.75 

June 21, 74 

65.35 

67.88 

85.99 

Jan. 5, 75 

71.02 

73.90 

93.94 
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Figure 4.5. Mobile Bay PCM (Oct. 1972) 
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Table 4.12. Joint Histogram Between GTM and Piecewise 
Linear Classification Map 
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Table 4.13. Similarity Measures of PCM's w. r. t. GTM 


Date 

Similarity Measures (%) 

Normal 

Boundaries 

Ignored 

Inventory 

Oct. 17, 72 

60.45 

63.15 

81.30 

Nov. 17, 73 

70.08 

72.65 

87.37 

Dec. 5, 73 

73.42 

76.32 

97.34 

Apr. 10, 74 

70.91 

73.69 

87.55 

June 21, 74 

68.41 

71.00 

96.19 

Jan. 5, 75 

69.76 

72.52 

90.02 
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Figure 4.6. Mobile Bay SNEHACM (Nov. 1973) 

37 


Table 4. 14. Joint Histogram Between GTM and Nearest 
Neighbor Classification Map 
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Table 4.15. Similarity Measures of SNEHACM’s w. r. t, GTM 


Date 

Similarity Measures (%) 

Normal 

Boundaries 

Ignored 

Inventory 

Oct. 17, 72 

59.27 

61.97 

80.78 

Nov. 17, 73 

72.56 

75.51 

88.55 

Dec. 5, 73 

66.22 

69.00 

87.39 

April 10, 74 

55.28 

57.74 

73.36 

June 21, 74 

58.71 

60.65 

82.36 

Jan. 5, 75 

55.18 

57.31 

78.02 
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using estimates of means and covariance matrices for the various classes based 
on training samples. The result of this method is called PRITICM. Figure 4.7 
shows a PRITICM of the June 1974 data set, the joint histogram of which appears 
in Table 4.16. The similarity measures (w.r.t. GTM) of the PRITICM’s corre- 
sponding to all the six data sets are presented in Table 4. 17. 

4 . 5 Comparisons Between Classification Maps 

It can be seen from the joint histograms in Section 4.4 that, in many cases, 
the off-diagonal elements are larger than the corresponding diagonal elements, 
particularly in the case of classes other than forest (3) and water (4) . This is 
due to the fact that forest and water classes have large homogeneous regions 
while the other classes have relatively large numbers of boundary pixels, and 
errors in registration of the GTM and the Landsat images have the most signifi- 
cant effect near the boundary regions. To support the conclusion that the dis- 
similarities found between the GTM and the CM*s are attributable to registration 
errors, several of the CM's for the December 1973 data set were compared 
relative to each other. The corresponding similarity measures are shown in 
Table 4.18. These numbers should be compared with the "normal” columns in 
the previous tables . It can be seen that these are much larger than the similarity 
measures relative to the GTM. Also, the joint histograms (not shown here) are 
all found to have dominant diagonals. 

4. 6 Chi-Square Tests for a Typical Case 

Chi-Square tests were applied to test several hypotheses in the case of the 
comparison between the GTM and the LCM for the October 1972 data set. Starting 
with the joint histogram (contingency table) the values were computed as de- 
fined in [2] for nine different hypotheses listed below: 

1. The distribution of the classification inventory agrees with the 
distribution of the ground truth inventory. 

2. The distribution of the correctly classified pixels agrees with 
the distribution of the ground truth inventory. 

3. The distribution of the number of correctly and incorrectly 
classified pixels is optimum with respect to the given in- 
ventory and without regard to class. 

4. The correctly classified pixels are randomly distributed. 

5. Each ground truth feature is randomly distributed among 
the classification features according to the ground truth 
inventory . 
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Figure 4.7, Mobile Bay PRITICM (June 1974) 
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Table 4.16. Joint Histogram Between GTM and Histogram-based 
Maximum Likelihood Classification Map 
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Table 4.17. Similarity Measures of PIUTICM's w. r. t. GTM 


Date 

Similarity Measures (%) 

Normal 

Boundaries 

Ignored 

Inventory 

Oct. 17, 72 

62.98 

65.88 

88.68 

Nov. 17, 73 

71.22 

73.80 

91,49 

Dec. 5, 73 

73.45 

76.31 

93.68 

April 10, 74 

70.75 

73.60 

90.39 

June 21, 74 

70.77 1 

73.45 

90.44 

Jan. 5, 75 

69.01 

72.00 

92.82 
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Table 4.18. Similarity Measures between CM's for Dec. 73 

Data Set 
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6. Each classification feature is randomly distributed among the 
ground truth features according to the classification inventory. 

7 . The distribution of the number of correctly and incorrectly 
classified pixels is random without regard to class, 

8. The numbers of correctly and incorrectly classified pixels 
for a particular class are randomly distributed. 

9. The distribution of the classified pixels is independent of 
the ground truth. 

The results of the test are shown in table 4,19. Here, the test numbers corre- 
spond to the hypothesis numbers above and D. O. F. is the number of degrees 
of freedom. It is found that the values are all large and the corresponding 
probabilities are very close to zero indicating that the above hypotheses are 
false. However, these tests are too stringent [2] . 

4.7 General Remarks 

4.7.1 Qualifiers 

The following qualifiers are necessary in drawing the conclusions from 
the data presented in this chapter. 

• The GTM used here was a preliminary draft accompanied by a 
disclaimer stating that it was meant only for field checking and 
review. However, it was used in the comparisons here and 
this step should therefore be considered preliminary. 

« The UTM coordinates of the control points needed for the slight 
correction in orientation of the GTM were obtained from a map 
of scale 1:125000 with a similar disclaimer. 

• The Landsat data were corrected using affine transformations 
to UTM coordinates, the parameters being obtained based on 
several ground control points. While the GCP's were chosen 
sufficiently scattered over the scene and the rms errors were 
smaller than one pixel at the GCP's, the peak errors were of 
the order of 2 pixels at the GCP's and could be 3 to 4 pixels 
elsewhere. 

4.7.2 Conclusions 

f It was found that the point by point C'normal") similarity measures 
w. r, t, GTM in the case of the data sets considered here are in the 
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Table 4.19. Results of the Tests 
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same range as those for the Bald Knob, Tennessee data set [1] . The 
results for HINDUCM, LCM and ETCM are slightly better than in [1] . 
These are the only techniques common between the present ease and 
[1] and they support the intuitive conclusion that large homogeneous 
regions tend to increase classification accuracy. 

• The data sets from the November and December 1973 scenes result 
in higher similarity measures than the others in the case of most of 
the techniques considered. Also, the October and June scenes seem 
to result in lower similarity measures in all cases . It is found from 
the classification maps that these scenes have more cloud cover than 
the others and, since the cloud pixels are not identified and removed 
from consideration before finding the similarity measures, they con- 
tribute to misclassifications'. 

• The similarity measures between classification maps of the December 
data set are significantly higher than those between the CM’s and GTM. 
This could be due to the fact that the CM's are perfectly registered 
relative to each other while the registration errors contribute to 
disparities between the CM's and GTM. 

• The percentage of boundary pixels is smaller in the present data 
sets compared to that in [l] and the increase in similarity measures 
when the boundary pixels are ignored is not as significant. But, 
when registration is improved, it is expected that there will be 
greater increases in these siinilariiy measures. 

t In [l] it was seen that the similarity measures based on inventories 
only were very close to those when the boundary errors were ignored. 
It is found in the present case that the inventory similarities are much 
higher. This, again, points to registration errors. 

• The results of the histogram-based, supervised, table look-up tech- 
niques (SNEHACM, PRITICM and PCM) depend on the grid-size 
chosen in deriving the histogram. Here the same grid-sizes as 
obtained in the case of HINDUCM' s of the corresponding data sets 
were used for the supervised CM's. Whether higher similarity 
measures can be obtained with finer grids (i.e. higher radio metric 
resolution) needs further examination. From the point of view of 
computation time, however, it is not feasible to use these tech- 
niques with grid size unify (i.e. full resolution that is available in 
the data). 
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4.7.3 Suggestions for Further Work 

It is clear that an improvement in registration of the GTM and the Landsat 
data sets is needed before further comparisons are made. When the final version 
of the GTM is received, the changes, if any, should be incorporated and an up- 
dated digital version of the GTM produced. 

The classification maps in the present work have been at level I. It is 
difficult to obtain reliable training samples for a level n classification. There- 
fore, it seems unnecessary to have a GTM at level H. A significant reduction 
in the manual effort involved in establishing the correspondence between the 
region numbers and class numbers on the GTM would result if a level I GTM 
were available. The present digital GTM is at level n and, since much of the 
manual work has already been done, it seems appropriate to finish further 
corrections of it at level It. These corrections should include modifying the 
locations presently having -l*s (confusion due to missing boundaries between 
distinct ground truth classes) and 99*s (indicating that labels were unavailable 
on the GTM supplied) . 

The GCP (ground control point) coordinates should be determined more 
accurately and a nonlinear (instead of the affine) transformation be found for each 
of the data sets. It may not be sufficient to use just the EN term. Better results 
could be obtained (when the GCP's are accurately located) with all the quadratic 
terms included in- the transformation. It is desirable that the error in registration 
be less than one pixel throughout the image. 

The reasoning behind obtaining similarity, measures with the boundary 
pixels ignored is that the boundary pixels tend to be mixtures of classes and 
hence errors in classifying them should be "excused". Thus, only the errors 
in classifying the Interior pixels of a homogeneous region are attributable to the 
defficiencies of a classification technique. This argument can be extended to the 
case where there are uncertainties in registration. If the maximum error is 
known to be N pixels then the subset of the image to be considered for comparison 
should consist only of the points farther than N pixels from the nearest inter- 
class boundary. H such a subset is large enough to be statistically significant, 
it will prove to be useful in comparing the classification techniques. 
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5. 


TEMPOHAL CHANGE DETECTION 


5.1 Introduction 

The repetitive, eighteen day cycle in the coverage of each Lands at satellite 
was planned in order to have the capability of examining changes in the same scene 
as a function of time. Hence it is advisable to examine some concepts and metho- 
dology for defining and detecting temporal changes through processing of multiple 
data sets collected in different Landsat passes over the same scene. It is assumed 
in this chapter that the data sets have previously been spatially registered. 


Temporal change detection as a concept is easily defined as the process of 
identifying the changes that have taken place in a scene. However, the method- 
ology to be adopted is to be tailored to meet the needs of the user in terms of 
what precisely is the change the user is looking for. For example, the user 
could be looking for a change in the shape of certain features in the scene, like 
a body of water, or the emphasis may be just on the inventory of certain classes, 
say, forest acreage. Thus, different view points and needs of users have to be 
accommodated in an automatic information processing system. It is therefore 
desirable to first establish the identify of each individual pixel in all the data - 
sets uniquely which may then be further processed to delineate the type and 
extent of change the user is interested in. Further, the correspondence between 
the clusters or classes of the different data sets has to be established through 
inventory /spectral matching or manually through supervision. 

Once the pixel identify in each image is established and a correspondence 
between the labels in the two images is derived, the problem is then reduced 
to overla 5 dng the two images and flagging the pixels appropriately to identify 
the nature of change occuring at each pixel as well as whether the pixel is in 
the interior or boundary of a region. The details of the processing involved 
and the methodology developed to cater to this need are presented in the following 
sections. 

5.2 Overview of the System 

The information flow in the system, schematically represented in Figure 
5.1, can be briefly described as follows. With the two data sets representing the 
two temporally separated images of the same scene as input, the preprocessor, 
consisting of a classification/ clustering package (depending on whether the environ- 
ment is supervised or unsupervised) and a registration and geometric correction 
package, produces two classification maps GCMl and 2 which are geometrically 
corrected and registered with each other. The first image is then processed by 
the software subsystem FLGBDIES to identify and flag the interior and boundaries 
of clusters, leading to what will hereafter be referred to as the reference map. 
Either manually or through automated inventory / spectral (e.g. MXSMLKTY) 
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Difference Image Producing Subsystem 



Figure 5.1. Schematic Representation of Temporal Change Detection Methodology 
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matching schemes, a cluster or class correspondence list correlating to two 
classification maps is derived. This list along with the reference map and the 
second map is then input to the MASKER subsystem. Here, a difference image 
of the two input images is produced which will have in coded form all the' infor- 
mation (concerning each pixel) necessary for recognition of all possible types of 
changes occuring in the scene. In addition, the coded difference image has also 
the necessary information to determine whether a given pixel is an interior or 
boundary pixel in the reference map. With this coded image as input, the CHange 
Recognition And DEpiction (CHARADE) Systeniis now ready to process the user's 
input request for detection and depiction of any particular type of change in the 
scene, as for example, change in the boundary of a certain class, say, the water 
bodies in the scene, etc. 

5 . 3 Description of the Change Detection System 

As described earlier in the overview of the system, and portrayed in Figure 
5.1, the change detection system consists of three major components (subsystems) 

• A preprocessing subsystem 

e A difference image producing subsystem 

• A change recognition and depiction subsystem 

The preprocessing subsystem in turn consists of several smaller packages . 
The first significant constituent of the preprocessing subsystem is a classification/ 
clustering scheme (depending, of course, on whether the environment is super- 
vised or unsupervised) which, given the two temporally distinguished data sets 
(images) produces two classification (cluster) maps. The classification tool 
could be one of the several evaluated and reported elsewhere. These classifi- 
cation maps are then registered either relative to a common reference such 
as the UTM coordinates or with respect to each other. This is done by identify- 
ing several ground control points and finding appropriate geometric transfor- 
mations (using the program GEOGREF) . The transformations are applied to 
the two images using the program GEOCOR and the resulting images are shifted 
relative to each other to ensure proper overlaying. 

The registered maps are then compared by using either automated inven- 
tory matching scheme such as MXSMLRTY or spectral matching or through the 
available supervision in the environment to derive a class (cluster) correspon- 
dence list correlating the class (cluster) number in the two maps. Thus, the 
output of this preprocessing system mil be the two geometrically corrected 
maps GCMl and GCM2 and a class correspondence list LCC . 
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The difference image producing subsystem consists of two components: 
a boundary flagging scheme which identifies and flags the boundary pixels (a 
boundary pixel being defined as a pixel with at least one neighbor different 
from itself) and a difference image coding scheme which produces a difference 
image coded so as to label uniquely the type of change occuring at the pixel. 

In addition, the code also identifies the pixel either as interior or boundary 
pixel. This code can be written as 

IPIXEL (•) =1IX (•) -1)*M+LCC (IY(. )) 

where IX (•) is a value ranging from 1 to M for interior pixels and M+2 to 2M+1 
for boundary pixels, depending on the class to which the pixel is assigned in 
Map 1 with M as the number of classes in Map 1, lY (•) is a value varying from 
1 to N depending on the class to which the pixel is assigned in Map 2 with N 
as the number of classes in Map 2, and LCC (•) is the class correspondence list 
with numbers ranging from 1 to M designating the equivalence between the classes 
in Map -2 and Map 1 . 

Thus the coded difference image represents a completely processed infor- 
mation encoded image which can be employed to determine and depict the type 
of change of interest to the user. 

The CHAnge Recognition And DEpiction (CHARADE) subsystem essentially 
consists of a process designed to cater to the user's requests for depiction of a 
particular type of change. The user's input request code is decoded and the coded 
difference image is scanned to identify and categorize the pixels belonging to the 
class(es) of interest according to whether they are 

• Pixels undergoing the type of change of interest to the user 

• Pixels undergoing no change 

• Pixels undergoing changes, but not the type specified by the user. 

In addition, the rest of the scene of interest is flagged to distinguish it from the 
areas outside the image. The inventory of the different categories as well as 
a coded map depicting the pixels belonging to these different categories are pro- 
duced and recorded. The output images arc written out on tape and, in addition, 
printer plots are also displayed. 

Thus, the change detection system, given the two raultispectral data sets 
corresponding to two temporally distinguished images of a single scene, leads 
to change depiction maps of interest to the user. 
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5.4 Experimental Results 


The change detection system was implemented on IBM 360/65 and tested 
using the multispectral data sets corresponding to the three passes of Landsat 
over Mobile Bay area on December 5, 1973; June 21, 1974 and January 5, 1975, 
The change detection scheme was applied to identify the changes from December 
to June and December to January. The classification maps were produced by the 
linear sequential method. Depiction of six different types of change was requested: 

• total changes in all classes. 

• changes in only the interior regions of all classes. 

• changes in only the boundaries of all classes. 

• additions (or gross growth) to class 1 (identified manually 
as Urban. 

e deletions (or gross reduction) from the interior region of class 3 
(identified manually as Forest) 

• changes in the boundaries of class 4 (identified manually as 
water) . 

The statistics of these changes are given in Table 5, 1 

Change depiction maps are shown in Figures 5.2 -5.6 for the December- 
January pair of classification maps. 


5.5 Discussion of Results and Concluding Remarlcs 

As can be observed in Table 5.1, more of the changes occur along the 
boundaries than in the interior regions , This is only to be expected as boundary 
pixels are not only more susceptible to changes but also likely to be mixtures 
of more than one land use class or category. The figures for inter-seasonal 
and intra-seasonal changes have to be compared and interpreted in the light of 
the fact the intra-seasonal case deals with winter season data sets only, al- 
though involving a temporal separation of 13 months. On the other hand, the 
inter-seasonal case while involving a temporal separaticn of 6 months only, deals 
with winter and summer data. The statistics reported in Table 5.1 confirm that 
the inter-seasonal changes are far more significant than other temporal effects . 
This is especially true in the case of Forest, for example, where the change from 
December 1973 to June 1974 is significantly hi^er than the change from December 
1973 to January 1975 although the temporal separation is a lot more. Though only a 
few particular types of changes were looked into in these experiments the system 
has considerably more flexibility in depicting various other types of changes as 
discussed earlier. 
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.Table 5 , 1 . ,Temporal..C.hangeJDetectionJRe6ults 



Type of 
CH^ge 

DATA SETS CONSIDERED 

No. 

Decembers, 1973' 
INTER-SE 

Vs, June 21, 1974 
ASONAL 

■ '■■■ — _ - . 

December 5, 1973 Vs. Jan. 5, 1975 
INTRA-SEASONAL 



No . of pixels changed 

Percentage change 

No. of pixels changed 

Percentage change 

1. 

Total changes 
in all classes 

442013 

31.99 

322484 

23.24 

2. 

Changes in 
interior regions 
and ail classes 

158100 

18.22 

73566 

8,42 

3. 

Changes in the 
boundaries of 
all classes 

283913 

55.24 

248918 

48.44 

4. 

' Additions to 
class 1 
(URBAN) 

80745 

5. 84 

67516 

4.87 

5. 

Deletions from 
the interior of 
class 3 
(FOREST) 

112692 

33.49 

i 

j 

37738 1 

- - 1 

11.22 

6. 

Changes in the 
boundaries of 
class 4 
(WATER) 

4500 

40.77 

• - — 1 

4938 

1 

^ 

44.60 
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Figure 5,2. Total changes in all classes (represented by dark pixels) 
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Figure 5.3. Total changes in interior regions (black = change, gray=no change, 
white =boundary) 
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Figure 5.4. Total changes in boundary regions (black = change, gray = no change, 
white = interior) 
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Figure 5.5. Additions to the urban class (black = change to urban, g^ay = no change, 
white = change to other than urban) 
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Figure 5.6. Deletions from the interior of the forest class (black -deletion 
from forest, gray = no change, white =pixels other than interior 
of forest) 
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DATA COMPRESSION 


6.1 Introduction 

Remote sensing by means of multispectral scanners such as those employed 
on the Landsat program and those planned for the Earth Observatory Satellite can 
result in extremely large quantities of data. The Landsat scanner, for example, 
could generate approximately 125 reels of 1600 bpi magnetic tape per day. It is 
apparent that data compression would yield benefits in the recording, transmission 
and storage of this information. Some of the more obvious benefits are reduced 
on-board storage, simpler data transmission, reduced ground data recording, 
and fewer data tapes to archive. 

Data compression is accomplished by exploiting the structure or redundancy 
which exists between data samples. This means that the data does not take on 
all values with equal probability or that the value at any point is not totally in- 
dependent of the data at every other point. It is apparent that spatial correlation 
exists, due to the extension of generally well defined regions on the ground over 
several neighboring pixels. Hence, the purpose of data compression is to trans- 
form the image data in order to reduce the degree of correlation between the 
samples so that redundancy in transmission is minimized. 

6.2 The Karhunen-Loeve Transform 

The Karhunen-Loeve [l, 2] transform (also known as the eigenvector trans- 
formation, the principal components transformation, or the Hotelling [ 3 ] trans- 
formation) results in uncorrelated transform samples, thus minimizing redun- 
dancy. This transformation is optimal, if the criterion used is the mean squared 
error between the original image and the reconstructed image. The matrix re- 
quired for implementing the KL transform has as its i*ows the eigenvectors of the 
covariance matrix of the data. It then follows that the principal components are 
uncorrelated (the covariance matrix of the principal components is diagonal), 
and the variances are the eigenvalues of the original covariance matrix. 

Disadvantages are that the correlation matrix of the image must be known, 
necessitating two passes through the data, and a time-consuming matrix multi- 
plication must be performed. 

In order to describe mathematically the transform process, an image is 
modelled as a sequence of discrete data values with certain statistical properties. 
Experimental evidence indicates that a large variety of imagery data can be re- 
presented by a so-called markov sequence of variables having correlation de- 
creasing exponentially with distance between terms. The autocorrelation matrix 
R is then given by 

Rij=P^^"^L 0 <p<l. 
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^N-1^ ^ samples, the elements of R range in value from unity to 

p . The correlation between neighboring samples is shown in the following 
definition: 


x.i = PXi_i + £i_i, 

where the sequence {e^} has zero mean and variance l-p 2 . The KL transform of 
the sequence is defined as the transformation that diagonalizes R, with 
elements the eigenvectors of R. Thus the transformed sequence {xf} has un- 
correlated samples with variances given by the eigenvalues. 

The KL transform of the sequence [kj] is given by [4] 


U X 2 1 - p2 

where Xj = 1- — and 

1 - 2p cos + p2 

[o)j j are the positive roots of 

tan NU, = . 

cos 0) - 2p + p2 cos u) 


Since the values of u) are so defined (by a transcendental equation), the transfor- 
mation is defined by nonperiodic sine functions, and no computational simplifi- 
cations are possible. 

6.3 A Fast Karhunen-Loeve Transformation 

The KL transformation may be obtained from the Fast Fourier Transform 
(FFT) if a proper modification of the data is applied. [5] This is the minimum 
variance or interpolative representation, which is defined as 


X,i=Xi+ Vf 

where is a linear combination ofjhe remaining terms in the sequence, i.e. , 

{Xj, jj^i } . If the variance, (xj ~^i » minimized, the definition of x: 

reduces to ^ 

Vi)+ vj. 
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Since this definition involves terms Xq and the boundary -or end conditions 

must be defined: 

Xq = PXi + Vq 
XN+i = PXN + vn+1 

The minimum variance, sequence does have nearest-neighbor correlation, and 
hence its correlation matrix (NXN, tridagonal) has the form B^Q where 



The minimum variance representation may be written in the following 
manner: 
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The coefficients are the elements of the rows of matrix Q, and so the equations 
may be written in matrix form as 

QX = V + B 

where B is an NX I vector containing only the information at the end points, b, 
and bfj, with all other elements zero. Thus there is established a correspondence 
between the original sequence {xj} and the minimum variance sequence {vj} . 
furthermore, the eigenvectors of Q are [6] 





N+1 


sin 


ij.TT 

N+1 


and the eigenvalues are 



cos 


3 tt 

N+1 ‘ 


The W form the KL transform matrix of the sequence [Vj3 (being the eigenvectors 
of Q, the correlation matrix), and are periodic sine terms, and thus computable 
by an FFT algorithm. 


Applying W to the matrix equation for {xj^} > we obtain 


or 


WQX = WV + WB, 


A 

X 



using WQ = X. W and denoting the KL transform WX by X, etc. For zero mean 
data, the boundary conditions may be approximated by zero, and the KL trans- 
form becomes 



Extension to two-dimensional image data is readily obtained, fe ] The repre- 
sentation 


which defines xj in terms of its two nearest neighbors, is extended to include its •• 
eight nearest neighbors, as shown in Figure 6. 1. 
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Figure 6.1. Coefficients of eight Nearest Neighbors, 
where a = — • 

l+p2 

6.4 Implementation of the Algorithm 

The steps required for coding a data compression program via the KL trans- 
form are as follows; C?] 

(i) Calculate the statistics (i.e. mean, variance, horizontal and 
vertical correlation parameters) of the source image. 

(ii) Create a zero mean image, by subtracting the image mean 
(or mean of the image block in the case of block by block 
coding) from each point in tiae image. 

(iii) Use the equation involving the data sample and its eight 
nearest neighbors to compute Vi j. 
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(iv) Take the two dimensional- sine transforms of vy to obtain vj j . 

^ /V 

(v) Calculate xy = v.j / ^ y • 

A 

(vi) Quantize xy using n^j bits to obtain the values to be used for 
transmission. The number of quantization levels required for 
transmission is proportional to the transmitted value, xy , and 
hence proportional to l/Xy (vy being relatively small and 
constant) . The number of bits required to transmit 1/ \y is: 

log2 (1 Ay) = - log2 ^ y ♦ 

Hence, the number of bits to be assigned to the component (ij) varies according 
to 


njj = bi - b2 logg \y . 


In terms of the. number of bits assigned to the major component, and m, 
the average number of bits/pixel, ny is given by: 




= nii- 


(nil - m) logHj/ij^j^ 

2 Z log\y/\ij^ 
i=l 3=1 


Since this expression must be converted to an integer for each nj[j, the actual 
bit rate obtained varies somewhat from the predicted value, m. 


6.5 Results 

Two computer programs described in reference [7 ] were implemented on a 
255 X 255 segment of Landsat data covering the city of Mobile. The programs 
were an image analysis program and a Fast KL transform coding program. As 
supplied, the IBM 360/65 CPU times required were 6 l/2 minutes and 3 1/2 min- 
utes, respectively. 

The image analysis program performs the following tasks: 

(i ) Compute and print the histogram and statistical parameters 
of the image 

(li) Calculate the horizontal and vertical correlation parameters 
over all image blocks . 
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(iii) Compute and list a desired bit rate constant (m) versus actual bit 
rate table . 

The Fast KL transform coding program does the following jobs: 

(i) Create differential image {vy] of a 15 x 15 image block. 

(ii) Apply Fast KL transform to {vy} . 

(iii) Calculate bit assignments to different elements in the transform 
domain. 

(iv) Perform quantization. 

(v) Apply inverse Fast KL transform. 

(vi) Store final result as a 255 x 255 image on magnetic tape. 

(vii) Compute and print Iiistogram and statistics of tlie encoded image. 

The total run time of 10 minutes for a 255 x 255 image is long, but is readily 
reduced by adjusting the quantization of the desired bit rate table, and by saving 
the image analysis parameters required in the transform coding program. 

The statistical parameters of the original image and the reconstructed image 
are given in Table 6.1. The histograms of the input and final images and of the pixel 
by -pixel differences are given in Figures 6.2, 6. *3, and 6.4. The original and re- 
constructed images are shown in Figure 6.5. 
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Table 6.1. Statistics of the Input and Output Images 



Input Image 

Output Image 

Minimum 

13,00 

10.07 

Maximum 

82.00 

67.01 

Range 

69.00 

56.94 

Mean 

20.63 

20.50 

Standard 



Deviation 

4.00 

i 

i 

4.03 

» 

Mean Squared Error Between Images 2 . 34 
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Figure 6.2. Histogram of Input Image 
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Figure 6,3. HistogTam of Output Image 
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Figure 6.4. Histogram of Output Image Minus Input Image 



KARHUNEN-LOEUE TRANSFORM 



ORIGINAL IMAGE 



1 BIT/PIXEL RECONSTRUCTED IMAGE 


Figure 6.5. Original and Reconstructured Images of the City of Mobile. 
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PAKT n 


1. INTRODUCTION 

This part of the report Is a formal documentation of the programs devel- 
oped for the analysis and evaluation of multivariate decision methods for classi- 
fication of remotely sensed data and change detection. 

There are ten sections in this part, each section documenting one major 
software element. The programs are not detailed at the subroutine level, but 
ar e 'explaihed in terms of the inputs and outputs that one needs to know as a user. 
The subroutines needed for satisfying the external references are tabulated in 
each case. 

The programs and most of the subroutines are in FORTRAN IV and are 
implemented on an IBM 360 with the H compiler. They are all available as load 
modules on a users’ library. The names of the data sets on which the programs 
documented here were located at the time this report was prepared are; 

SMART .DASARATY . LIBRARY 

SMART .RAMPRIYA'. D091576 . LIBRARY 
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2. GEOMETRIC COEHECTION 

2.1 NAME 
GEOCOR8 

2.2 PURPOSE 

To apply geometric correction using nearest neighbor rule to a large rec- 
tangular image. The transformation from the output to input coordinate system 
can have eight parameters, six of them accounting for rotation, scale change, 
and shift and two providing a second degree term with the product of the output 
coordinates. 

2 . 3 CALLING SEQUENCE 

This is a main program. It is on a partitioned data set as a load module. 
The member name is GEOCOR8 . 

2.4 INPUT -OUTPUT 
2.4.1 Input 

The following input parameters should be supplied in data cards according 
to the formats and read statements indicated below. 

READ 100, NREC, NEL 

READ 200, A, XO, YO, ALFA, BETA, SX, SY 
100 FORMAT (216) 

200 FORMAT (6F12.3) 

where 

NREC, NEL are the number of records and the number of pixels per record 
in the input image; 

A is a 2x2 matrix accounting for rotation, scale change and skew; 

XO, YO are the shift parameters; 

ALFA, BETA are coefficients of the product term; 

SX, SY are scale factors. 
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The transformation applied would then be 


'x' 

- Ap 

n 

T" 

-xo-j ^ 

'ALFA' 

_Y_ 

[yi 

> 

_YoJ 

BETA 


XP*YP 


XP~ 


XPP/SX" 

YP_ 


_YPP/SY_ 


where 


XPP, YPP are the coordinates in the output image and 
X, Y are the coordinates in the input image. 

Note: Generally, this program is used with A, XO, YO, ALFA, BETA 
found using a mean squared error minimization process with ground control 
points (e. g. GEOGREF [1]) . The units of A are input pixels per km. Then 
SX and SY will indicate the number of output pixels desired per km in the XP 
and YP directions respectively. 

The input image data should be as unformatted FORTRAN records with 
NEL words per record and one pixel per word. 

2.4.2 Output 

The output image data will have NRECO records (unformatted FORTRAN) 
with NELO words per record and one pixel per word. The values of NRECO and 
NELO are printed along with the coordinates of the top-left and bottom-right 
corners of the image, (see Section 2.9, Method). Also, some details about the 
implementation of the program are printed. 

2.4.3 File Storage 

This program requires a direct access file of 1500x1500 bytes for ii iterme - 
^iate storage ( NREC s: 1500, NELs 1500; for larger values of iMeC or NEL, 
the DEFINE FILE statement and the space allocation for unit 90 should be changed) . 

2.5 EXITS 

Not applicable 
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2.6 USAG E " 


The program is in FORTKAN IV and implemented on the IBM 360 using the 
H compiler. The program is in the users’ library as a load module. 


'2.7 EXTERNAL INTERFACES 


This program calls several routines . The linkage is indicated in the 
following table . 


Calling Program 

Programs Called 

GEOCOR8 

SARN 

DAWN 

GEOMM 

GEOMM 

i 

GEOMl 

DARN 

VNATS 

GEOM2 

GEOM3 

SAWN 

1 

GEOMl 

VMOV 

GEOMIB 

XCHNGE 

GEOM3 

1 

GEOM4 

GEOMIB 

GEOMIA 

GEOM4 

DARN 


2 . 8 Performance Specifications 
2.8.1 Storage 

This program is 192680 b 3 d:es long, mainly due to an array DC di me nsioned 
48000 words . This array can be reduced in size, but the cost is an increase in 
direct access reading. Including the external references and buffers this program 
needs 256 K bytes of storage. 
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2.8.2 Execution Time 


The time required depends largely on the output image size which in turn 
depends on the input image size and the transformation parameters . The time 
needed to correct a 1200 x 1200 Landsat image to UTM coordinates to produce 20 
pixels per km, thus generating approximately 2150x1850' pixels of ou tput, is about 
14 minutes on an IBM 360/65. 

2.8.3 I/O Load 

None except as specified by Section 2.4. 

2.8.4 Hestrictions 

NREC ^ 1500; NEL S1500 (See Section 2.4.3). The numbers on the input 
image should be between 0 and 255. 

2.9 Method 

The details of the method are presented in Section 2.5 of part I. These 
steps are implemented by the routine GEOMM. The main program GEOCOR8 first 
reads the input data from a sequential data set (unit 10), converts them into bytes 
and copies to the direct access data set (unit 90) . The processed image will 
appear as a sequential data set on unit 8. The main feature of GEOMM is that 
it requires only one work array IX which it allocates for various buffers de- 
pending on the computed Val ue o f the output record length. The details of the 
subroutines follow the description in section 2.5 and are also apparent in the 
comments in the attached listing. 

2.10 Comments 

This program is designed to handle resampling with the nearest nei^ibor 
rule but can be modified easily to perform bilinear or bicubic interpolation. 

The present method of data handling involves considerably more I/O than that 
described in [3] which used segmentation, but was designed because computations 
needed for the nonlinear transformation would be complex under that approach. 

2.11 Listings 

The listings of this program and the associated routines are attached at 
the end of this section. 
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reproducibility of the 
ORIGINAL PAGE POOR 



2.12 Tests 


The program has been cheeked out using a test pattern, applying a 45° 
rotation to it and printing the results . Also, a transformation with a small non- 
linear term added to the 45° rotation has been tried. The program has been 
used to correct a classification map of the Landsat data of the Mobile Bay, Ala. 
test site to UTM coordinates using a nonlinear transformation. 
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LEVCL 21*8 ( JUN,JL^JL 


05/,36C.., FORTRAN H 


OATF 76,2S?/20>r5449 


I5N QQ02 
ISN 0003 
ISN 033^ 
ISN 0005 
ASU ODQii_ 
ISN 0007 
ISN oooa. 


AHE.5 MAJllii a £.T^02>UN£ CN.IeS6.i.5.LI£5QCCHK4 

SOURCF,RflCDlCfNQLISf|NnDECK,LOADfHAP,NaeD!T,IDfNnXREF - 

_CaMMC]N/GEaH/A,XQ#yQi ALFA, BETA ,J)£T1 ,DET2 , 1LD,IH1.,JLQ, JHI . B ,NRE AD. 

DIMENSION AC2, 2) ,0(2,21 

-DIMENSION „IX(ABOGC) > 

EQUIVALENCEUXI AOnn,LX( in 

I nG!CAL»i-.u(J-i3 cinx 

DATA HAXC/49''0^/ 

-DEF1NE-.FILE 90 ( ISOD, 1 50C, L, IA3tl 


ISN 0009 


-READ. INPUT. IMAGE..J] 
READ KO,NRFC,NEL ' 


ISN 0313- 
ISN 0011 
ISN Q012 
ISN 0013 
J-SN,-aaiA, 
ISN 0^15 


COPY INPUT IMAGE TG DISK* 

-NELAS.NELPA 

DO 10 I^lfNREC 

X ALL -5ARN.I IC,1X,NELA) 

DG 20 J=1,NEL 

gjujirixi j) 

CALL 0AWN(9P,I,LX,NEL 1 


READ TRANSFORMATION PARAMETERS. 

J HE-.E I GHT ..EAR A H ,.X0 ^ YU , ALFJV , 0 ETA_ARE .THQS L-EOU_ND_ 0 Y...G,EO.GR££^_ 

SX AND SY ARE SCALE FACTORS IN THE E AND N DIRECTIONS IN PIXELS/KM. 


^ ISN 0017 
. ISN -3318-.- 
ISN 0019 

^lShUQ32Q--. 

ISN 0321 
lSN..QDZ2-_ 


ISN 0023 
-LSN-2a2^-.. 
;SN 0025 
JSN. .0026 
ISN 0027 


-^C MQ D.I F-Y: JTH E.. XRANi^GRldA UDN_A.CC_0UNT 1 NG JEOR^S CALE-JEJAITDRS.* 

A(l,l ) = A(), ll/SX 

AU,1J=A( 2aJ./5X..- 

An,21^A( 1,21/SY 

A I_2 , Zl=JVi.2,.21/-S.Y 

ALFA=ALFA/SX/SY 

: B.E-UHi£XA25X/J5.y 

r 

_C AaPLY-jrHE-TRANSXQRKAtlQK.. l : 

CALL GEOMM( IX,MAXC,NREC,NEL,8I 

SXO£ 

100 FQRHAT(2I61 

-Xinn EQRKAii^iaa^ai : 

END 


fepBODPCl 


LEVEL 21.3. ( JUN-.I4_J_ 


0$/360 PQOTRAN H 


ISN 0002 
ISN 0003 
.15 N. 030.4. 
ISN 0005 


>TinNS “ H 4fiE* HAIN. gPT-C Z ,L1N£CN.I =56^ 5 J Z.EiOCOOK.t 

SOURCE. EBCDIC, NOLI ST, NODE CK .LOAD, HAP, NOeO I T» ID, NQXREF 

.SUBROUTINE GEOHH 1 1 y ,H AXC , NR EC ,NEL , NTAPU I . . 

DIMENSION IXIMAXCT 

j; OMHQ N/ii E OH/ A , XQ ,Jf 0 , ALFA,.B £.TA,JltIl , DE T2 , 1 LO ,i HIjJLQ, JH I . B . NEE A D_ 
DIMENSION A(2,2) ,6(2,3) 



purpose: TO APPLY GEOMETRIC TRANSFORMATION TO A LARGE IMAGE STORED 

DN..A..DIRECT.. ACCESS DATA..SET(UN1,T_9C). .THE OUTPU.T. .HI LL APPEAR ON 

SEQUENTIAL ACCESS DEVICE NTAPO. 

aNL)L.DNE ARRAY. IX..QF. MAXC H(IR0S_IS SUPPLIEC AS_A, WORK. AREI A.N.0. .THji_ 
ALLOCATIONS FOR VARIOUS SUBARRAVS IS HANOLFO INTERNALLY. 


C inputs: a, XU, YO, ALFA, BETA 

.£ IJLA.NS£ORMALiaN_OJE_T.HE_F.f>RM s_. 

C (X Y) = (XP YP)A» ♦ (XO 

J: WHERE .THE DENOTES TRANSPOSl 

C INPUT IMAGE, XP,YP ARE THOSE 

X THF IMAGE ON UNIT 9r » 

C 

_C £1ULPUIS‘!.._.A ,.XO,...YD,_ALFA, BE 

C FOR CONVENIENCE OF TMPLEMENTA 

J: ia_B_SQ_T H A J_LH E;_S I VE N _.T R A N SF 

ILO, IHI, JLO, JHI ARE THE 

!ja AC£ ).„0 F „TH E _.TOP_ IND . 0 OH 

C THE OUTPUT IMAGE. 

C the_output„image hill ..have.. 

C WORDS PER RFC'ORDIQNE WORD 


ARE EIGHT PARAMETERS DEFINING A 

~Y0)~ ♦'(Alfa" BE'fAlTx'mp) 

TION, X,Y ARE THE.COOROINATES IN THE. 
IN THE OUTPUT IMAGE. 

I 


T.LARE modified. B1_THE„R0UTIN' GEOM.l 
TION. OUT, FIRST, A IS TRANSFERRED 
GRMATICN CAN PE _RECOyERED. IF.NEEOJ.D.i 
COORDINATES! IN THE TRANSFORMFO 
CH_RCH S_AND LEF T AND RIGHT COLUMNS 0 

_IHI-ILn*l RECORDS WITH JHI-JLO*! 

PER PIXEL). 


ISN 0006 
. ISN OQ07 


DEFINE FILE 9^ I NPEC.NEL ,L , I AV ) 


FIRST, COMPUTE THF OUTPUT IMAGE SUE AND MODIFY THE" SUPPLIED TRANS- 

_£ ORHALION_f OR . C O.N V E N 1 ENC E 0 F_ Ui£L £ HEN T A T I ON . 

CALL GEOH11NREC,NEL,NRECO,NELO) 

_HBJT£(6,10'?JNReC..NEL,NREC.a,MLC ; ; 

IN COMPUTING AN OUTPU T RECC R O. WE N FE O TO ALL OW FO.R._2,_«_MLP_K’O RnS_ 
FOR COMPUTING COORDINATES IN THE INPUT IMAGE CORRESPONOI NC TO EACH 
.POINT OF OUTPUT, NELO WORDS FOR. STORING THE OUTPUT. VALUES AND NELO, 
WORDS FOR A "CODE ARRAY”. THUS, HAXC-4*NEL0 WORDS ARE AVAILAbLE 

FOR A CI RCULAR B UFF ER FOR TH E INPU T RECORDS A ND A BUF FER POINTER 

ARRAY. ” 


ISN 0003 
ISN 0739 
ISN 0710 


MAXCP=HAXC-NFL0*4 

NR=HAXC£.?4/( N EL.+.4 1 

HRTTE(6,2rO)MAXC,KAXCP,NR 


FIND STARTING ADDRESSES FOB WORK AREAS. 


ISN 0711 


X COORDINATES 
IA01 = HAXCP+1 .. 


Y C OORDINATES. 

IAD2=IA01+NEL0 


ISN 0012 



LEVEL 21.8_t JUN._16_1. 


QS/36C . FOXTRAN H 


76.253/?'>.r6.<T^ 


ISN 0002 . 
ISN 0703 
..ISN_QQQA_ 


LL1QNS_=_MHE= B » Q£J^2 .U.N.EO t=f 6.. S I ZE^COO'^K . . . 

S^URCF, EPCOICfNOLI ST, NrPEC K ,L Cftn , HAP , NQFDI T , I f) ,N.3XRFF 
.SUBROUTtME CEOMllNREC,NEL,M’ECC,NELnj 

COHHON/GECM/A,XO,YO, ALFA,RETA ,DET1 ,0ET2 ,ILO,IHI , JLO, JHIte.NOEAO 

.DJM£NSIQN_.A<2,2I ,B(2,.2J. 

THIS ROUTINE MODIFIES THF GIVEN TRANSFORMATION FOR CONVENIENCE OF 

jCQMEUr.AXiaM,i , 

IT IS ARRANGED SUCH THAT THE NUMBER OF INPUT RECORDS TO BE HELD IN 
_CQRE_F.DR.-CDMPUTlNG.OHt-nUTPUT_ ftECDRb‘.IS_MlNlMLZEIL.JINO...J(, .Y INCREASi 
WITH RESPECT TO XP, VP RESPECTIVELY. 


ISN 0005 
ISN 0J06_... 
ISN 0007 

ISN 0003 

ISN 0007 

-ISN joau 

ISN 0011 


PRINT THE INPUT PARAMETERS. 


HRITEI6, no) 

_H Rl.TE.Ii> ,2CC M I A I L, J 1 , J?.l , 2 J ,I I 

WR!TE(6,3''OIXO,YO,ALFA,0ETA 

_CALL_VMOV(A,A,BI L 

0ET1»AC1 , 11 ‘■AI 2 ,2)-A(l,?» ‘■AI 2,1» 

-D.EJ.2s. A U UJ.tDE.LAr.A I 2, H»ALEA 

CALL GE0M1B(NREC,NELI 


ISN 0012 
ISN 0013 
ISN OOIA . 
ISN 0015 
ISN 0016.. 


FIND MAXIMUM NUMBER OF INPUT RECORDS NEEDED TO PRODUCE A RECORD OF 
.0UTP.UTLR1...II1.THE. GIVEN TRANSFORMATION IS USED! J12_IF...TH£ ,R0LES_0£_ 
XP AND YP ARE INTERCHANGED). 

.Rl.rABSUA(l,2)*ALFARIL0)p( JHL=lJLO) I 

Rl»AHAXl(Rl,ABS I(A( l.aUALFAPIHnPIJHI-JLGU) 

.A2 = ABS(IAtl,lUALFA«>JLCI)*nHl-IL0H 1 

R2=AHAXU R2,ABSl(An, l)*ALFA*JHI)*nHI-rLtn) 

JIJ11BI.LE.R2JGQ _T.0.JIC 


ISN 0018 

ISN 0019 

ISN 0020 
ISN 002L . 
ISN 0022 


ISN 0023 
ISN 0025 
ISN 0026 
.ISN.-312i. 
ISN 0023 
ISN 0029. 
ISN 0030 
ISN 0031 
ISN 0032 


CALL XCHNGE(A(l,U,A( 1,2)) 

_CALL_XCHNGE1AI2,U,AI 2,2») 

CALL XCHNGE(ILO,JLG) 

.£ALL..XCHNCE(IHI.» Jhl) 

CONTINUE 

"h odT FY The transformation, if necessary, toi ensure that 

-WITH XP, , 

IFU( 1, 1 ) .CE.''. )G0 to 2P 

-AU,n=rAIl,l) .. . 

A(2,l )=-A( 2,1) 

-4 LEA.? -.ALFA 

0ETA=-OETA 

-CALL- XCIr!NCE.U.LQjJ_H-l) . 

ILO^-ItO 

JH!=-1HI_ 

CONTINUE 


INCREASES 


ISN 0033. 
ISN 0035. 
ISN 0036 


MODIFY THF TRANSFORMATION, IF NECESSARY, TO ENSURE THAT V INCREASES 

H1TH_VP. _ _ 

IF(M2, 2).GE.b. )G0 TO 30 

..A.U»2»E_-A(1,2) .. 

A12,2)«-A(2,2I 


C 

C 



PA6I: n>'2 


ft L FAa rJLEA . ! 

BETA=-ftETA 

CALL. XCHNCEtJLO»JHn __ . 

JLQ=-JLO 

JHIa-JHI i •_ 

30 CONTINUE 

NR.ECQs.lHl^ll.QtJ 

NELO SJH1-JLO+-1 

~c PR I N Tmo d I f I e6 ’ tr an sf or h a t i'o?j 

MRITEI6, AC9) ..J! 

HRITE(6,20n){(A(I,JI,J=l,2l,Ial,2l 

URITP(A.3 0''tVn.Vn.AI CA.HRTA 


WRITE(6,5POIILO,JLO,IHI,JHI 

RETURN , 

100 FORMAT!/// • GIVEN TRANSFORHATioN PARAtii TE RSI M 

-200 F0RHATI//-'_MATRIX'/(1X2E15.7H 

300 FORMAT!/' SHIFT VECTOR : • 2E 1 5. 7/ 

t rnpci: ir T CMTC nc oanpiirT tcdu^iipic ii 


400 FORMAT!/./' MODIFIED TRANSFORMATION PAR AMETERS :• I 
-SftO EUfillAll//d_IQR.. LEFT,.CCR.NEB?Lllti*.'I6i'). .. BOITOH RIGHI.C! 


ISN..D0A5. 




uni 2 U 8 .1 JUHJTA^L 


0 S/ 36 C ^ FORTRAN H 


OATP 7 S. 2 '^ 3 / 2 '^.r 6 . 10 


ISN 0002 
ISN 0)03 
l$U COO^. 


mDNi^N 3 ^ME=_MAIN#QPT» 02 .LINECNT = 56 . 5 IZE= 0 C:i:K*i 

SOURCE, FOCniCtNOLISTfNOnECK ,l QAO , MAP, NOEOI T, 1 .0 ,NnxRFF 
..SUORQUTINE GFDMlAUtY,XP,YP,ISIGN) . . . 

COMMON/GFOM/A,Xa,Ya,ALFA,BETA,DETl ,OET2 ,I LO ,! HI , JLO, JHI ;b ,NRE AO 
^DIMENSION A12, 21 ,6(2,21 , . . . .. 


tSN 0005 
ISN CDOb. 
ISN 0007 
ISN OOOd. 
ISN 0009 


ISN aon 

ISN 0013. 
ISN OTIA 
ISN .0015 


■JJNIL.-X2, YR.JU-VZ.N.X,.,.Y, 

XXO=X-XO 

BB=DETUALFA<‘YYO-BETA«'Xxb 

.CCrYYD^^Al 1.2 )tXXQ?A(2,2) 

I F( A6S(0ET2)*GT. l.E-aiXP^I-QB^ISIGN^-SQHTI BB^fiB-4«OET2«'CO 1/ 
l?onFT 2 l 


IF( ABS(DET2).LE •1.5-31 XP=-CC/BB 

Jt P^(.Y-Y.O- &.L2 >.Uj» tA» XP I 

RETURN 



LEVEL 2l.8„(, JUt4 7^ ) 


CS/16^ FpftTPAM 


date ,7ft.2S3/2?.f6..1A„ 


ISN 0002. 
ISN 0103 
ISN 0304- 


TIONS - ( iJAME= ti,ft,lNtOPT;l?»UNEC<)IT,=56|?UP-^0O'>Kt 

S OURCE I f PCD K, NPU I ST, NOOECK, LOAD, HA P, NOE niT,!r),‘JQX«EF 
.SUPRrUTINF CEOH10(NREC,NELI 

COHHCN/CEDH/A,xa,Yr', ALFA.OETA ,0E T1 ,0E T2 , t LO , I H I , JLO, JH I , a , NRE AD 

JDI.HE.NSION.. A.(2,2) ,B(2,21 


_E.1N CL-M I N_ .4 H!1_KA X_VJiJL U E S DP X P . VP F CB., .X , .Y „RA NCltlC_JlR OH L.TQ_NREC. 
AND I TO NFL RF SPECT t VELY. 

_C»LL_GEQH!A«1.,1..,XP.YP,1> 

CALL GEOKlAd. , 1. ,XPH,YP!1,-1I 

..l.SlGN'l L- _• 

IF(ABS(XPH>.6E.4BStXPHC0T0 10 

_L5J£li£iJ 

XP*XPH 

~xm'in=xp ' 

-.xHAx^xp : 

YHIN=YP 

Y H ft Ks V P 

NREci^HAXCfNREC-l ,1 J 

_NELl?HAXC(NEL-l,l) _ 

DO 20 K«l,2 

~iEL=NErr' 

_JF(K.EQ..llGQ.JJQ_aH 

IREC=NRECl 

"continue 

_P.Q_iQ_J=JLi N P EC.,.IJ5J£ 

00 20 J»1,NEL,1EL 

_X=.l 

Y = J 

-CALL GEOMlAlX,Y,XP.,.YP-,.l.SI.Gfil _ . 

XMAX=AHAXUXMAX ,XPI 

_XHIN=AMIN1(XH!N,XP». ; 

YHAX=AHAX1(YHAX,YP) 

_YHIW=A HTN1( YHIN..YP1 

CONTINUE 

- ILOrXHIN-1,5 ... ... 

IHl=XHAX+1.5 

-J[-LQ=lYMIN-.1..5 

JHI=YMAX«-1, 5 




CDHRILER-DETIGHS - tJAHEs_.JiaiN.DP-Lsl2.«i,lMfCNT==5fi ,iI2E ■=CRHaJU 

SOURCE, EBCDIC ,NOUSr,NOOFCK,LllflO,MflP,NOEf)IT,ID,NOXREP 

.. I5N 0002 SUBROUTIME CEOM 2 1 X ,Y . I , ME L C,NREC , ^EL^ 

ISN 0003 C0HM0N/CECM/«,X0,Y0,4LFft,BETA,CETt ,OET2,IIO,IH1 ,JL0,JHl,8,MREaci 

„ ISN 0004 OIKENSIQN A I 2, 2 » ,B( 2, 2» • 

C 

C F.lHD.-..ARRAVS-nF X ANO Y mnftOt NATES POB T HF t«TH fUlTPIIT R Pf gp.Q^ 

ISN 0005 OIHENSION XiNELOl ,Y(NELO» 

ISN 0006 11=I*.ILD::1. 

ISN 0007 00 lO. J»1,NEL0 

ISNOOOfl JJ = J+JL0-1 J 

ISN 0009 X(J)=A(l, 1,2I*JJ+XC*AIFA6II»JJ 

IsK-OOiii 

ISN 0311 IFIl.GT.XUl.OR.XIJI.CT.NRECICa TO 10 

ISN 0013 Y.UI3.AU, l.L6XLAA(2,2UJJ+Y0*aETA*Il»JJ„ 

ISN OOU IFn.GT.VUl.OR.YI J>.GT.NELIYU»«r, 

ISN 0016 10 CONLINUE 

ISN 0017 RETURN 



. LEVEL 21.8 ( JUN .7^..L 


CS/36C PC»TFAN H 


. Q M" f . - . 76. g 51/20.»'6.?A 


jOmpjlEILa PJ..LIMS-^ £*_.a A IN . Q PI « 0 2 , LIUECM T » £6 , S I Z E *: CC.?. t ». 


SnilBCE, FPCOir rNPL ! SI.NOOFrK.l f1 AO , MAP , ^OEOI T, in ,NP XBFF 

ISN aQ02_. _• SUDROUTINE CECP. I (IN , X , Y , IREmc , I Ci ICB, Lf. , N.;F C , NFL ,t;FLO , IR I . I « F I 

C 

C COMPUTE, QNE RECORD OF R.ESAMPLEQ.OUTPUT OY READING. THE NECESSARY 

C • INPUT BFCGROS FROM THE PIRFCT ACCESS OEVICF (UNIT 9*' I . 

C IH15-RRDGRAH„JJ_EESJGNEQ_TD_HjaNDl.E NEAREST HFl f.HRn R AND BILIMFAr’ 

C- INTERPOLATION FOR RFSAMPLING. THE COMPUTATIONS OF JRI, JRF SHOULD 

. C DE._CHANCED IN. THE CASE OF BICUBIC INTERPOLATION.,— ' . 

ISN COO 3 COHHON/GEO'1/A,XO,YQ, ALFA, BETA ,DET1 ,DFT? ,IL0,1HI ,JLO, JHI ,B ,NRE AD 

— ISN 0004 DIMENSION. AI2,2»,B(2, 21 

ISN CD05 LOCICALfl IN 

— „1 SN n a o . 6 (UHENSION TCn(NR>,TREHC(NELOI ,X( WELPI tYINELP 1 . 1 Nf NFL .N R I . 101 N ELQ) 

C DEFINE FILE .90 I NRFC,NEL,L , lAV I 

C INI.T.1.4UZ.E IREHC. 

ISN 0007 DO 10 1=1, NELQ 

__ ISN.03Q8 IREHCni.= I . _.. 

ISN cao9 loi i»-*r 

LSK_CDJJ3 LQ LELY 1 II . E Q. 0 . 1 1 RFHf M 1 

C 

C J= INCLHAXl.-AND_Hl!!i..CF„BECORD.NUHai&.S..NEEDE,Q..TO .CD.HP.UTJE_XHE .REMAI.NIM 

• C PART OF THE OUTPUT RECORD. 

ISN 0012 IPASSeO 

ISN 0013 -30 JRIslorOf'OO 

ISN. 00 14 JBF.*a 

ISN 0015 IPASS=IPASS4l 

_ - ISN 0016 D0.2C .L=1.,NELQ 

00 ISN 0017 IFdREMCm.NE.ltGO TO 

_Oi _ C 

C IREMCIIl.FO.l INDICATES THAT I0( H ' HA s'nOT B 6 eF WMPUT Eo" Y E 77 

I iN-DOXa IKl5J( UJ 

ISN 0020 JR.I=MINr( JRItlXl » 

.. ISN 0021 JRF.=MAXO.URF,lxn.. . . ' •, 

ISN 0022 20 CONTINUE 

..ISN 0023 JEE3HINC(JRF41,NREC) • ; 

ISN 0024 IFIJRI.GT.JRFIRETURN 

C 

C NOW, JRI, JflF ARE THE HIN AND MAX RECORD NUMBERS N^EOFO FOR 

... C. XOMPUTING THE REMAINING PART OF THE OUTPUT RECORD. 

C IRI THROUGH IRF ARE THE RECORPS IN CORE. EXCEPT IN THE CASE 

. __C LP.ASStI,. as MUCH„AS POSSIBLE OF THF OUTPUT RFCORO HILL HAVF BFEN 

C COMPUTED USING IRI THROUGH IRF BEFCRE EACH CALL TO CE0M4. 

c 

C IFIIRI.LE.JRI.AND.IRF.GF.JRFJ’ THE FNTtPF OUTPUT RECORD CAN BE 

- C COH.P.UTED. 

ISN 0026 101R=0 

C Tf~T^Tl T . j’rF read forward' ANO~COHFUTE PART OF* THE (jUTPUT RECORd."” 

,.- 1S.>1J3.02.7 IFHRF.LT .JRFIl DlRd 

C 

C IF_I,RF,CE.JflF.ANO.IRI.GT.JRI JEA0.6ACKKAR0 AND COMPUTE PART OF 

C THE OUTPUT RECORD. 

ISN 00?9 ^IFtlR£.GE.JRF,AND.IRI.r.T,JRinCIR--l 

C 


PAGE 


ISN 0131 
ISN 0033 

ISN 0036. 

ISN 0037 


C IF...I£ASS^.£Q • LIH£fi»JJ.ai_TiiE-LlRSX. CALL _TD C EDH4,CnftRES P OMP INC TQ 

C the present record. THER'^FQREf COMPUTE AS MUCH AS POSSIBLE WITH- 

C OUT. READING ANY NEW RECORDS. . 

IFnPASS.EQ.niOlR^^ 

^ CaIl GE0M4( IO!R,NLFFTrIRl t IRF^NIRrJRIf JRF#lN|NEL7lCBiiRfMC,NFL0t 

. ia.x» y 1 

IFCNLEFT.CT.OIgO to 3C 

REJ.URN 

END 



Lf:Vet JUN..TA--L 


OS/36^ FOkTR^N H 


hate 76*2F3/2^>^6>?9 


ISN„0:372 

. ISN 0003, 

ISN 030^ 

.J,5M„Q0a5 

ISN 0106 


Jl£IJD,NS^r^N^Me?- MAlN.0PT = 02iUNECNT^56 , 5 IZ E:^"OC:!K» 

SOURCE , FdCpiC ,N0LIST,NGDf CKtLUAOrf'APfNOFO!Ti rO^UJxnF 

5JH^R0UTIN£ GEO'^Al lOIR ,NLE F T 1 1 R1 , 1 RF ,NR t JR I f JR F, IN,NEL i jC A ^ I ^ f 

» NELOtTCfXfYI 

LOGICAL^i IN 

DIMENSION INiNELfNR), ICBCNRlflREHCCNFLOl I XlNELOlfY(NELC) , iniNELni 

COMHON/OEDH/A,xa,Va,ALPA^eETA tCETlTDEfF^TiaTlHI , JLO, J h7, bTnR E AD 

c1^17tTT’pART of of the DUlPuf RFCORO 8Y RFAOING ir^HANY RECORDS 

AS_PQSSIBLF..QF.THE„.INPUT^.IHAGE^IMHEOIATELY AHEAO^OF^OR »6HIN0 

THOSE ALREADY IN CORE* 

IH IJU£RaGKAM-Jii^Dm.GJ^EO^F.Q.R^£.AR£lLJ£lGHDDa...iNIlREilUU.O/lt_J:^ 


ISN OJQZ NB2JL5NEl$L2ri 

ISN OJOd IRIO^IRI 

ISN 0109. I FnDIR)lOi20,30 


ISN ono 

ISN 001l„. 
ISN 0C12 

ISN oon._.. 

ISN OIU 


.00 ISN 0015. 
00 ISN 0016 
ISN O017_ 
ISN 0018 

_^JSJy_001j9___ 


ISN ''lao 
ISN 0321 
ISN 0022 


ISN 0023 
ISN 0)2A 
ISN 0)2§.. 


ISN 1026 
ISN 0*^27, 
ISN 0023 
ISN (?030„ 
ISN 0031 
ISN 0033 
ISK 0036 
ISN 0035 
ISN 0036 
ISN 0037 
ISN 0033 


IRl=MAXOnRI-NR*l,JRU 
,IRF = MINOlIRIfNRTliJRF). 
IRl-IRI 

_iB2riRio-i : 

GO TO 60 


IDiR.GT.r* FIND IRl, IRF TO RFAD FCRWARD* 

..IRFQ«JRE^ : 

IRF*:HINr( IRF^NR-lf JRFI 

.JK.UMAXO( mF=jyRU.tJSl 1 ... 

IR1=IRF0*1 

.JJR2^ia.E 


_MQai£Y_ICD SUCH. THAT, KBCD.GmS- ADDRESS OF .mJLIR-XNPUT. 
INC*IRI-IR1C ♦NR21 

_DQ.50..I?lfNR L 

ICBIIUMODC rCBI n*INC, NR) + 1 


READ THE NECESSARY INPUT RFCOROS, 

„Dg 60 I^IRl,IR2 . . . .. 

NRFAD=NREAD+1 

_CALL_DARNI0r,i^, IN.Il.ICP(I--IRimi#NFL ) 

Q MP.U1E_P A R TS O F THE OU TPU T ARRAY T HAT C A N P E COMPUTED^ _ 

NLEFT=‘p' ' 

_Pn.J0.l£l*NFLO 

IFUREHCI n.NE.DGO TO 7 ^ 

_iETXlILt.5 . ^ , 

IF( IR.UT. IR I.QR . IR.GT.IRFKO TO 75 

TR^TC B( IR-TRUl ) : 

K=Y(n^.5 

__igmniNi iCriRi _ - . 

lRF.MCm = 2 

_ 6D^T0 70 ..... , .. 

NLEFT^NLFFT+l 


RECDRO., 



OFlT: 


MGE_C''? 


ISK 0339 1ft (LQHUNUE— 

IS^ ft9<tO RETURN 

• ISN 03-il ' END - . . 


LEVEL 21.9 t JUN „ CS/36C '‘O'- TRAM H .. OATF. . 76, 25.3/2r* ''±•3* 

C.PaeiLEIL-aE.T.lItMS. rt. NAM Es MAIN , aPI.s.0 2 ♦ LINECN T = ?6 , S U E =CC ?CR.. 

SQURCE, E0Cr)rC,NOLIST,NOOECK,LO4Pt‘'AP,NOFOIT, lO.'nvRPF ' 

ISN 0002 . SUBROUTINE VNATSCIX.N) ■ . . * 

tSN 0103 DIMENSION IXIN) 

ISN OIDA DO_.lQ_Ls.l,N 

tSN 0105 in ixm»i 

' ISN OJO& ELEIUaU 

ISN 0007 


END 




.LEVEL 21.8.1 JUN..7.AJ .CS/36C FOftTPAN H ' ... . ... .DATE. .7,t..252JL2.ajuCjfej.ia. 

CDrtP-LL£R-Q£XJCttS-.~-M.A.M E = MA.IN * D P I.=£.2L. LUtECtLL?- 56 ..5 IZ C CiJL# L„ - ! 

SCURCE,EPCDK.NnLIST,NaDEi:K,Lniin,HAP,NnEDIT,ID,NnXHeF 

ISN 0002 SUaRDUTINE XCHNGE(X,YJ ... 

ISN 0003 H=X 

ISN 0005 Y=H 

IS7J_flaa6 SLEIUfiil 

ISN 0007 


END 




3. THUSHNING OF BOUNDARY IMAGES- 

3.1 NAME 
PEELS 

3.2 PURPOSE 

Starting with, the output of a micro densitometer digitizing a boundary image, 
to apply a given threshold of density and reduce the thickness of the boundary lines 
by "peeling" their outer layers while preserving the distinctness of regions separated 
by them. 

3.3 CALLING SEQUENCE 

CALL PEELS (NTAPI, NTAPO, NREC, NEL, IT, MPASS, MDEV, 

NDEV, LX, LY, IBDY) 


where 


NTAPI, NTAPO are the logical unit numbers of the input and output sequen- 
tial data sets; 

NBEC, NEL are the number of records and the number of pixels (bytes) 
per record in the input image; 


IT is a threshold on density; if IT is positive (negative) all points with 
densities 2: IT (^ IT ) will be regarded as boundary points; 

MPASS is the maximum number of iterations permitted (see Section 3.9, 
Method) ; 

MDEV, NDEV are logical unit numbers of two direct access scratch 
data sets defined as indicated in the listing of PEELS; 

LX, LY, IBDY are scratch arrays with LX, LY dimensioned as indicated 
in the listing and IBDY dimensioned NEL . 

3,4 INPUr-OUTPUT 

3.4.1 Input 

The input image should be on a sequential data set with unit number 
NTAPI and consist of NBEC records and NEL bytes per record, each record 
corresponding to a line of the digitized image and each byte, to a pixel. All other 
inputs are as indicated in the calling sequence. 
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3.4.2 Output 


The output of this program will be on unit NTAPO as a sequential data 
set with NHEC records . The records will be in SLIC (scan line intersection code) 
format. That is, the first word of the I'th record indicates the number of words 
that follow and each subsequent word is a column coordinate of the intersection 
of the I'th scan line with the boundary image. 

3.4.3 File Storage 

This program requires two direct access scratch data sets to handle the 
intermediate iterations of the boundary data. The sizes of these data sets are 
indicated in the listings attached. 

3.5 EXITS 

No nonstandard exits. 

3 . 6 USAGE 

The program is in FORTRAN IV and implemented on the IBM 360 with the 
H compiler. The program is in the user’s library as a load module, 

3 . 7 EXTERNAL INTERFAC ES 

This subroutine calls several subroutines and the linkage is shown in the 
following table . 


3 . 8 PERFORMANCE SPECIFICATIONS 

3.8.1 Storage 

The subroutine PEELS is 1458 bytes long. However, including a driver 
(whose size depends largely on the dimensions of LX, LY, IBDY which are 
functions of NEL), &e re quired subroutines and the buffers the program needs 
approximately 70K for handling NEL =2100. 

3.8.2 Execution Time 

The execution time is highly dependent on the size and complexity of the 
boundary image, the thickness of the boundary lines and the maximum number 
of passes QVIPASS) requested. In the case of the Mobile Bay GTM (a 4000x2100 
level n map with boundaries 3 and 4 pixels thick) the initial thresholding, and 
reformatting took about 10 minutes and the subsequent interations about 6 min- 
utes each, with a final reformatting and cop3ung step taking about 7 minutes. Thus, 
with MPASS=4, it takes about 40 minutes of CPU time to process the image. 


93 



Calling Program 

Program Called 

PEELS 

PET 

SARN* 

VLTHR 

CMPRES 

DAWN* 

PEELER 

DARN 

EXPBDY 

CMPEES 

ISTORE+ 

. PEELER : 

i 

SVSCI 

PEELRl 

PEELRO 

DAWN* 

EXPBDY 1 

ILOAD+ 

i 

PEELEl 1 

i 

DARN 

BLSFTV 

BRSFTV® 

PEELRO 



icOMPl+ 

IAND+ 

BLSFTV 

— BLSFTV 

ILOAD"^ 

ISTORE+ 

BRSFTV 

ILOAD+ 

ISTORE+ 


* Entry under DAKN 

+ Logical function available in the user's 

library under a main member name LOGFUNC 

• Entry under BLSFTV 


94 






















3.8.3 Restrictions 


None 

3.9 METHOD 

The program has three major steps: 

(i) Thresholding, compressing and writing on a direct access unit. 

(11) Iterating to "peelV boundaries . 

(iii) Changing to SLIC format and writing on output sequential data set. 

3.9.1 Thresholding and Compressing 

The routine SARN reads each record (of NEL bytes) of the input data set 
into the array LX. The routine VLTHR thresholds each of the NEL bytes in LX. 
A logical vector LY is defined as follows: 

IF (IT.GE.O)LY(I) =LX(I).GEJT 
IT (IT.LT.O)LX(I).= LX(I).LE.IABS(IT) 

for 1 = 1, NEL. 


The routine CMPRES is then used to pack the information in LY into the 
first NEL bits of the array LX. The I'th bit of LX is ’’set” if and only if LY (I) 
is .TRUE. ' 


The compressed boundary information is then written on the direct access 
unit MDEV using the routine DAWN. 

3. 9 .-2 Iterating to Peel 

The main peeling routine is called PEELER. The input to this routine is 
from MDEV whenever IPASS, the iteration number, is odd and the output then will 
be written on NDEV. When IPASS is even, the inpiU; and output designations are 
interchanged. One call to PEELER removes one "layer” of the thick boundaries 
from top, left, bottom and right. 

To decide whether a particular boundary point should be deleted (i.e. the 
bit corresponding to it changed to 0), we examine a 3x3 neighborhood centered 
around the point. Consider the array 

a b c 
d e f 
g h i 
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where each letter represents a binary pixel. It is to be decided whether e, which 
is presently equal to 1 should be changed to 0. The conditions for a 'top peel' will 
be derived below and those for peeling from the othe^r directions follow by s 5 rrometry . 

First of all, e should be a top boundary point. That is, there should be 
no boundary point directly above e and there should be a boundary point below e. 
Therefore b = 0 and h=l are necessary conditions. Suppose bh=l, (Here, b de- 
notes the complement of b) . Then, we need only check whether e is a nonessential 
boundary point, that is, whether two O’s in the 3x3 array which are disconnected 
will stay disconnected where e.is made 0. Connectivity, in this context, is defined 
as the existence of a path not including I's and consisting only of horizontal and 
vertical segments. 

Now, it is ea^ to see that e is essential if and only if ad = 1 or cf=l. 
Therefore, the condition for a top peel is that 

bh (a+d) (c+f)=l. 

Equivalently, to perform a top peel we set 
e=e (b+h+ad+cf). 

It is co nven ient to'-f^lement the above equation by employing bit manipulation 
routines operating on pairs of 32 bit words, thereby performing the top-peel 
operation in parallel on 32 pixels . This is done by using the "current" array 
in place of e, the "previous" array for b, the "next" array in place of h._ Also, 
the previous, current, and next arrays are right (left) shifted by one bit .and- 
used for a, d and g (c, f and i) respectively in the peeling formulas . 

The routine PEELER minimizes the movement of data in core by using 
circular buffers for storing the "previous, current and next" arrays. An array 
J dimensioned 3 is used to store the indices pointing to these arrays (J(l) — 

previous, J(2) «»• current, J( 3) ~ — s^-next) and after finishing each record, 

only the array J is updated. 

Also, top, left, bottom and right peels are performed one after the other 
by just one pass through the data (thus minimizing I/O) by storing the intermediate 
results in core and operating with a phase lag. 

When the I'th record LX is read from the input data set (see PEELRi), 
BLSFTV and BRSFTV are used to generate arrays LXL and LXR with the bits 
in IiX shifted by one bit to the left and right, respectively. Next, the (I-l)th record 
is peeled from the top. The top-peeled output of the (I-2)nd record is peeled from 
the left. The top-and left-peeled output of the (I-3)rd record is peeled from the 
bottom. The top-, left- and bottom- peeled output of the (I-4)th record is right- 
peeled and written on the output data set. Also, whenever any peeling is done other 
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than from the right the output is shifted to the left and right by one bit and the 
results are stored in the appropriate core locations pointed by J(3), K+ 1. 

The routine PEELRO wth the appropriate ISIDE will perform the peeling 
of one record. The above operations performed for 1 = 1, NREC +4 will complete 
one iteration of peeling, constituting one call to PEELER. The number, NP, of 
words of input that were changed is counted during each call to PEELER. If NP = 0 
or the number of calls to PEELER has been MPASS, the iterations are stopped, 

3.9.3 Converting to SUC 

Each record is read from the last scratch unit on which the output image 
was created. The routine EXPBDY is used to sense each bit in the record. The 
bit number of each 1-bit is stored in IBDY. The total number, N, of 1-bits 
followed by N words of the array IBDY are written on unit NTAPO . 

3.10 COMMENTS 

On large images this program takes a long time to execute. To avoid loss 
of data on long runs it is suggested that the direct access data sets be saved 
so that, with slight modifications, the routine PEELS can continue where the last 
run stopped due to insufficient CPU time. 

3.11 LISTINGS 

The listings of PEELS and most of the associated routines are attached at 
the end of this section. The routines not included are: PET, a routine used for 
printing time elapsed between sections of a program; SVSCI, a routine which sets 
all elements of an array to a given constant; DARN and the associated entry points 
for a37ray read/write and the logical functions under member name LOGFUNC. 

3.12 TESTS 

The program was tested on a small portion of a boundary image, the image 
printed before and after peeling and was found to work satisfactorfly. 
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Vt-L I JAN n I 


OS/360 FORTRAN H 


DATE 76.ua/19.11.32 


COHP ILE R OPT IONS - NAHE« HA I N,GPT»02 , L INECNT «56i S IZ E=OOOOK, 

S OURC E * E BC D I C , NOLI S T . NODE CK , LO AO , HA P NOToTt, I 0,‘nJ XH E> 

tSN 0002 SUBROUTINE PEEL SINTI « NT O.NR EC |N£L, IT , HP AS S.HOE V»NOEVf LX .LV. tbOY ) 

ISN 0003 DIMENSION LXI U ,LY(l I , IbDYI 1 » 

C 

C DIMENSION LX(36>»nNEL>l)/B2^nitLY(|NEL-l)/44>li 
C DEFINE FILE HDE V (NREC » ( NEL-1 1/32+ 1 *U , I AW 1 1 
C DEFINE FILE NDE VI NREC« ( NEL>1 1/32+1 lUtl AV2) 

C 

ISN OOOA N»(NEL-1) /32+1 

ISN 0005 CALL PET(2I ’ • • 

ISN 0006 ■ DO '10 1*1 rNREC 

ISN 0007 . CALL SARNINTI ,LX,NEL) _ ' _ 

ISN 0008 CALL VLTHRUX ,NEL » 1 T, LY I “ ‘ 

ISN 0009 , , . LX(N>=0 

ISN 0010 CALL CMPRES(LY,NEL,LXI 

ISN roll 10 CALL OAWNIMOEV. l»LX,N*Al 

ISN C012 CALL PETI2I 

ISN DOIl . „ ,_D0 20 IPASS»1 iHPASS 

ISN COIA IF(MQD(IPASS,2I.EQ.1KALL PEELER(HOEVtNOEV»NREC*N»Li(Ux(li*N + l) » 

. LXIZA'+N+H .LYtNPI 

ISN 0016 • IFIMODI IPASS.21.E0,0»CALL PEELeR(N 0 eV,M 0 ev.NREC,M,LX.LXU 2 *N + ll, 

. LXI2A+N+II ,LY,NPI 

ISN 0018 • PRINT 100*IPASS,NP 

ISN 0019 CALL PETI2I _ _ _ 

ISN 0020 IFINP.EQ.OIOO TO 30 - — • - -■ 

ISN 0022 20 .... CONTINUE 

ISN C023 IPAS5*HPASS 

ISN 002A 30 JDBV*NDEV 

ISN 0025 "" IFIMODI IPASS, 2 l.eQ. 0 »J 0 eV-H 0 EV"' 

ISN 0027 DO AO 1=1 ,NREC 

ISN 0028 CALL DARN (JOEV, I ,LX ,N*A I 

ISN 0029 CALL E XPB OY I LX rN iNEL, I BDY» J I . . 

ISN 0030 AO WRITEINTOIJ,(IBOY(L)tL=l,JI 

ISN 0031 CALL PETI2I 

ISN 0032 ■ RETURN 

ISN C033 lOO FORMAT! 5X 'DURING PASS NUMBER'ISt* THROUGH PEELERJ!16»* WORDS OF COM 
-PRESSED BOUNDARY INFORMATION WERE CHANGED.*) 

ISN 003A END 



LEVFL 


JAN 73 


OS/360 FOATAAN H 


76,114/19.10.48 


COMPILER OPTIONS - NAME* MArN,0PT«02*LINECNT»56,S UE»OOOOK, 

SOURCE, EBCDIC, NOLIST, NQOECK.LOAOiMAP.NOEO'IT.TdVnOXREF ' 

rSN 0f>02 SUBROUTINE BLSFTVI IX ,N, I Y» 

!SN 0003 DIMENSION mN),lV(N) 

ISN QCD4 _N1»N-1 J 

ISN 0005 00 10 1 = 1 , N1 ■ ■ 

ISN 0006 lYI r) = ILOAO(IX( I<-1I,32,1» ■ 

ISN 0007 '10 ■ IVm = IST0RE(lX(n,im),32,3ir 

ISN 0008 IY(N»=0 

ISN 0C09 lYI N)=IST0RE(IXIN1,1Y(NI ,32 ,3l> . 

ISN COlO RETURN _ 

ISN 0011 ENTRY BRSFTVnX,N,IY» 

ISN 0012 _ lYm = ILOAO(IXl II ,32,31 I 

ISN 0013 00 20 1=2, N 

ISN 0014 _ IY(II»ILQA0UXm,32,31l _ _ 

ISN 0015 ' 20 lYl II»ISTQREIIKn-ll,IY(II,32,ll 

ISN 0016 RETURN 

ISN 0017 ' ■ END ' ' 
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LGVFL 21.7 ( JAN 73 } 


0S/3b0 FORTRAN H 


DATE 76.ua/19.10.57 


COMPILER OPTIONS - NAME* MA I N» 0PT»02 .LIMECNT *56 , S UE*0000k, 

SOURCE. EbCOICt NOLI STtNQOECKtLJ AO. MAP, NOEOIT* lO.NJXKEF 
ISN 0002 SUBROUTINE BPKINT (LX.N.NEL.LV.LO .LU ' 

ISN 0003 DIMENSION LXINI 

I SN OOOA LOGICAL I LOAD 

ISN 0005 LOGICAL*! LYINELI .LO.Ll 

ISN 0036 JHR0=1 

I SN CC37 JBIT=33 

ISN COOS DO 10 1*1 ,NEL 

ISN 0009 LY( n*LO 

ISN CCIO J6IT*J01T-1 

ISN OOU IFUBIT.NE.OIGO TO 20 

ISN 0013 J81T?32 

ISN 0014 JWRO=JHRb+l 

ISN 0015 20 IFdLOADlLXUURO) .JBIT.lULYm*!! 

ISN 0ri7 10 CONTINUE 

ISN 0018 PRINT 100, LY 

ISN 0019 100 F0RMAT11XA(32A1.1X)» 

ISN 0020 . .. RETURN 

ISN 0021 END < / 
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LbVEL 21.7 I JAN 73 1 OS/360 FOHTKAN H DATE 76.114/19.11.01 

CCmPILER OPTIONS - NAME- MAI N ,QPT*02 *L INECNT »56 »S Ue-OOOOK» 

■ SUORCE, EBCDIC, NOLIST, NODECK, L'JADrMAP.NOEO IT, ID, NJXREF “■ ’’ ' 


ISH 

CP02 


SUBROUTINE CMPR ES ( LX , NEL,LY 1 

ISN 

000 3 


LOGICAL*! LXINEH 

I5N 

0004 


DIMENSION LY(1) 

ISN 

CC05 


JHRO=] 

I SN 

0006 


JBITa33 

ISN 

0C07 


DO 10 I =1 ,NEL 

ISN 

0008 


J8I T»JO IT-1 

ISN 

0C09 


IF! jeiT.NE.OIGQ TO 20 

I SN 

0011 


JB1T=32 

ISN 

C012 


JMR0=JWRD*1 

ISN 

con 

20. . 

IX = LX(I 1 

ISN 

0014 


LYI JHRO) = ISTQREUX,LY( JUftDl ,JBIT,U 

I SN 

0015 

10 

CONTINUE 

ISN 

0016 


RETURN 

ISN 

0017 


END 




. REPRODUCIBILITY OP THE 
ORIGINAL PAGE IS IpOR 




LI-VFL 21.7 ( JAN 73 I 


OS/360 FORTRAN H 


DATE 76.1lA/19.11.0i> 


COMPILER OPTIONS - NAME* MA I N , OPT >02 , LINECNT *56 , SI ZE»CCOOK, 

SCJURCRf EBCOIC>NQU5TtNQDECAtLiJAU»HAP,hOEDIT»IOfNiXK£P 


I SN 

CC02 


SUBROUTINE EXPBOYUX»NtNEL,IBDY t 

ISH 

0003 


DIMENSION LXINI tlDOYm 

I SN 

0004 


LOGICAL ILOAO 

ISN 

000 5 


JWR0»1 

ISN 

0006 


JBIT*33 \ .... 

ISN 

000 7 


J = 0 

ISN 

0008 


. DO 10 I»l,NEL 

ISN 

0009 


JBIT«JDIT-1 

ISN 

mn 


IFUBIT.NE^OlGO TO 20 

ISN 

cci a 


JBIT=32 

ISN 

0013 


JWkO = jV)RD^l , 

ISN 

0014 

"■’20 

IFI.NOT.ILOADaXC JUROIfJBITairca TO 10 

ISN 

0016 


J*J+1 

ISN 

0017 


ieoY(j)*i 

ISN 

0018 

10 

CONTINUE 

ISN 

0019 


RETURN 

ISN 

0020 

. 

.. .END. . . . 
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LKVEL 21.7 C JAN 73 ) 


OS/360 FUKTKAN H 


DATt 76.n^/iy.lK08 


COMPILER OPTIONS 


ISN C002 
ISN 0003 
ISN onoA 
ISN OC05 
ISN C006 
ISN 0007 
ISN 0008 
ISN 0009 
ISN 0010 
ISN 0011 
ISN C012 
ISN 0013 
ISN 0014 
ISN 0015 
ISN 0G16 
ISN 0017 
ISN 0019 
ISN 0020 
ISN 0021 
ISN 0022 
ISN 0023 
ISN 0025 


ISN CC27 


ISN 0029 


ISN 0031 


ISN 0033 
ISN 0035 
ISN 0036 
ISN 0037 
ISN 0038 
ISN 0039 


• 20 


30 

10 


NAME* MAIN,UPT*02tUNeCNT«56,SIZE*OOOOK, 
SaURCtfEliCDICtNOLISTiNaDeCK.LOAOiHAPtNOEOITr lOiNJXREF 
SUBROUTINE PEELt R < M0£V t ND 6 V t NREC f N iLX t LXR t L aL# LY t NP I 
DIMENSION LX(Nt3t4) tLXR ( N t3 t A ULXU N « 3 |4 ) « L Y <N 1 1 J T3 » 
NR£C1*NREC+1 
NREC2«=NREC^2 
NREC3=NREC+3 
NREC4==NREC*4 
JdUl 
J(2>»2 
J(3)=3 

CALL SVSC ULXfN^'iatOI 

CALL SVSC I(LXR|12«J‘N,0I _ _ . 

CALL SVSCIUXL,12^N,0) 

NP^O 

00 10 I »l tNRECA 
DO 20 K«l,A 

• IFU.LE.NREC + KIGQ TO 20 

CALL SVSCKLXd t JOdKIfNrOI .. 

CALL SVSC niXR(l»JC3| fKUNiO) 

CALL SVScniXLI If J(3) tKI»N» 0 ) 

CONTINUE 

IFd.LE.NREOCALL PEELRl ( MDE V , I »LX , J ,N, LXRfLXL I 
IF( UGT.l.ANO.l.LE.NRECn 

• CALL PEElROCLXUflfnfLXRnflfllfLXLnflf llfJfNflf . , 

• LXil fJ(3)f2) fLXRUf Ji3l f2) fLXL (1 r J( 3 1 1 2 > 9 NPI 
IFI l♦GT.2.AND.I.Le.NReC2) 

« CALL PEEIROCLXI If 1 |2I 9LXRnfl92)fLXLC If lf2lfJfNf2t 

• LX(X fJ(3lf3IfLXR(ltJ(3lf3) 9 LXL ( 1 f J ( 3 1 1 3 1 1 NP I 
Ipn.GT.3.AND. I. LE.NREC31 

. CALL PEElKQ(LXn»l |31 fLXRIlflf3lfLXinf lf3l iJtNf3t 

• LX(1 fJ(3lf4) fLXRllf J(3) f4) fLXL U f J( 3 d A I f NP I 
IFII.GT.A) 

• CALL PEELRGI LXU 1 1 f A 1 1 LX R (1 f I f A » , LXU I « 1 # A I f Jf N fAf 

• LYfOfOfNP) , „ . : 

IFU.GT.AICALL OAWNl NDE Vf I-A f LY f A«N) 

DO 30 K«lf3 _ 

J(K)^MQO( JIKI f3l4>l 
CONTINUE 
RETURN 
END 
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LEVEL 2K7 ( JAN 73 I 


OS/36Q FORTRAN H 


DATE 76.11A/19.n.U 


COMPILER OPTIONS - NAME* MAI N t OP T-Q2 , LI NECNT *56 ^ S U E*OOOOK i 


ISM OC02 
ISN CC03 
ISM COCA 
ISN 0C05 
ISN 0006 
ISN 0007 
ISN 0008 


SOURCE,ettCDICtNOLIST,NOOECR,LOADtMAP,NOEDITtlDtNJXREF 
SUBROUTINE PE EL Rl ( NOE V , 1 , LX » J iNj LXR f LXL I 
DIMENSION LX(Nt3) t LXR ( N »3 I t LXL( Nt 3 1 1 J( 3 ) 

CALL OARN(NDEV»I ,LX(1,JC3II,N«4I 
CALL BL5FTV(LX(l,J<3lltNiLXLCL»;jnm 
CALL BRSFTVaX(l,jnn,N.LXRU,J(3m 
RETURN 
END 
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LEVEL 21.7 ( JAN 73 » 05/360 FflRT.RAN H DATE 76.114/19.11.22 


COMPILER OPTIONS - NAHE= MAIN.aPTs02,LINECNT=56,S IZE»ODOOK» 

SnURCe,EBCDtC,NULIST,NU0EC<V,L:)A0,HAP,N0E0IT,10,NUXKEF 
ISN C002 SUBROUTINE PEELROt LX , LXR ,LXL i J ,N, IS I OE t LY , LVR, LYL ,NP » 

ISN 0003 DIMENSION LX ( N.3 ) ,LXR ( N .3 1 . LXLI N. 3 I .LY ( N I , J ( 3 1 , lU ( 3> .LYRIN t »LYL IN » 

ISN P004 DO 60 I»1»N 

I SN 0005 LYI n»LXl I,J12» ) 

ISN 0C06 ‘ IFILYU ». EO.OIGD T.O 60 

ISN 0008 GO TO (10 »20i30t40> ,1 SIDE 

ISN 0009 10 Ium=>IQR(LXU.J(U) .ICOHPl (LXU.J(3l),32.32il 

ISN 0010 IH(2)sIAN0ILXRn«jm).ICOHPl(LXR(IiJ(2ll,32,32M 

ISN ecu NI 3 )cIaN 0(LXL(1 »jm)tlCaMPl(LXL(l.J(2)>t32«32)) 

ISN 0012 GO TO 50 

ISN 0013_ 20 „ U( U=IOR(LXR(I ,J(2I» .ICOHPIILXLU ,J(2n,32,32l>. 

ISN 0C14 IU(2)<:IAN0ILXR(I,jmi>lCQMPl(LX(I.J(lM.32,32n 

ISN 0015 lH(3)=lAN0ILXR(l»J(3n,IC0MPl(LX(l.J(3n»32«32n 

ISN 0016 GO TO 50 

ISN 0C17 30 Iwm<:I0R(LX(l|J(3n>ICaMPl(LX(I.J(ltl»J2,32l) 

ISN OCie IU(21<’IAN0ILXRn«J(3»)«ICUHPl(LXR(lfJI2n.32.32n 

ISN 0019 I«( 3» = lAN0(LXLU|Jnn.ICaHPl(LXL(l*JI21l,32,32H 

ISN 0020 CO TO eo 

ISN 0021 40 Iwm<IOK(LXLntJ(2l>»ICOMPl(LXRn,J(2»l.32f32l) 

ISN 0C22 ' (W(2)»IAND(LXLniJmi,ICQMPl(LXn.JUtl,32«32n 

ISN 0023 . Iri(3) = IANDILXL(I,JUM.ICUHPl(LX(I.J(3)).32,32}) 

ISN 0C24 50 iu( n-ioRiiHin tiwun 

ISN 0025 IU{ IMIQRUvHl) tlJOl I , . 

ISN 0026 LYI n-IANOILYin ilhllH 

ISN C027 IF{LXnfJ(2n.NE.LYm)NP»NP*l 

ISN 0029 60 CONTINUE 

ISN 0030 IF( ISIDE. EQ.4IRETURH 

ISN 0032 CALL 6LSF TV ( LY, N ,L YLI 

ISN 0033 CALL BRSF TV ( LY, N i LYRI 

ISN 0034 RETURN 

ISN 0035 . END 


LLVLL ^1.7 ( JAN 73 ) QA/360 FUKTKAN H DATE 76.114/19.11*38 


COMPILER OPTIONS - NAME* MAI N ,OPT «02 ♦ L INECNT *56 ,S I Z f •CiOCOK , 

SfH>RCF,EUCOICtNUUST.NilOeCK,LJA(),MAP,NUEOIT,IOtNIXREF 
ISN 0L>02 ■ SUBROUTINE VLTKR ( LX » N , I T ,LY I 

ISN 0003 LOGICALX'l LX( Nl tLY ( N) , F/. FALSE. /«T/. TRUE./ 

C 

C • THRESHOLO A VECTOR LX OF B BIT INTEGERS TO GET A T-F VECTOR. 

C IF IT.GE.n, LXdl.GE.IT IMPLIES LVm»T. IF IKO LXI 1 1 .LE. lAflSUTI, - 
C IMPLIES LYm»r. 

C 

1 SN 0004 ITT*IABSI ITI 

ISN 0005 IF( IT.LT.OIGO TO 10 

ISN rr07 DO 20 I»1 .N 

ISN 0008 _ LYm»F . . _ 

ISN 0009 20 IFILXn ).GE.lTT)LYm«T 

ISN roll RETURN 

ISN 0012 10 DO 30 I>1 >N 

ISN 0013 LYI I »*F 

ISN 0014 30 lFILXm.LE.ITTILYm*T 

ISN C016 . RETURN . . . 

ISN 0017 END 


O 

cx> 




FINDING DISCONTINUITIES IN BOUNDARY DATA 

4.1 NAME 
BOUDIM 

4.2 PURPOSE 

To find the discontinuities in digital curves stored in SLIC format. 

4 . 3 CALLING SEQUENCE 

CALL BOUDIM (IBDY, NTAPI, NREC, NEL, IRC, ND, NDIS) 


where 

IBDY is a scratch array to be dimensioned NEL *3 where NEL is the 
maximum number of boundary points in a given line; 

NTAPI is the logical unit number on which the input boundary data 
are stored; 

NREC is the number of lines (records) in the input data set; 

IRC is the output array of coordinates of the discontinuities; 

ND is the maximum number of discontinuities expected [IRC ' 
is dimensioned (ND, 2)]; 

NDIS is the output integer giving the actual number of dis- 
continuities found. 

NTAPI, NREC, NEL, ND are inputs to this routine IRC, NDIS are 
outputs . 

4.4 INPUT-OUTPUT 
4.4.1 Input 

The input data should be on logical unit NTAPI as a sequential data set 
consisting of NREC records . Each record should consist of the coordinates of 
the intersection of the corresponding scan line with the boundary image written 
as 


J, (X(I), I = 1,J) 

where J is the number of such intersections and IX (I) are the coordinates. 
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4.4.2 Output 


The output of this program is only through the calling sequence. 

4.4.3 File Storage 
None 

4.5 EXITS 

No nonstandard exits . 

4. 6 USAGE 

. This program is written in FORTEAN IV and is implemented on the IBM 360 
with the H compiler. It is available on the users* library as a load module. 

4.7 EXTERNAL INTERFACES 

The linkage with other subroutines needed with this routine is indicated in 
the following table. 


Calling Program 

Program(s) Called 

BOUDIM 

BOUDIS 

BOUDIS 

JCOUNT 


4.8 PERFORMANCE SPECIFICATIONS 

4.8.1 Storage 

This program is 834 bytes long. Including the external references listed 
above, the storage needed will be 2578 bytes (excluding the calling program which 
should provide storage for the arrays IRC and IBDY) . 

4.8.2 Execution Time 
TBD 

4.8.3 I/O Load 
None 
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4.8.4 Restrictions 


None 

4.9 METHOD 

The routine BOUDIM simply handles the I/O needed for finding the discon- 
tinuities. Connectivity, in this context, is defined in terms of the eight nearest 
neighbors of the point under consideration. Therefore, while examining the ith 
record of data, it is necessary to have the (i-l)st and (i+l)st records in core. 

The movement of data in core is avoided by using a circular buffer IBDY dimensioned 
(NEL, 3) and indexed by the pointer array IND dimensioned 3. Initially, IND is 
set to fi, 2, 3} . Always, IND (2) points to the current row. The numbers of 
boundary points in the three rows stored in core are in (NB(IND(J)), J = l, 3) . 

The routine BOUDIM starts by reading the first record into IBDY (*, 2) . Then, 
for 1 = 1, NREC the (I+l)st record is read into IBDY (*, IND(3)). The (NREC+l)th 
record is undefined. Therefore, in that case, NB (IND(3)) is simply set to 0. 

The routine BOUDIS is called to determine the coordinates of the discontinuities 
on the I’th record. Then the pointer array IND is updated. 

The functioning of BOUDIB is as follows . Each of the boundary points in 
the current record is treated as the point e in the following array. 

a b c 
d e f 
g h i 

The number of boundary points in this array excepting e is called the connectivity 
count of e . The connectivity count is calculated by examining the arrays IBDY 
f„ IND(2)), IBDY (*, IND(l)) and IBDY (*, IND(3)), stopping the calculations 
when the count equals 2. If the count is smaller than 2, then the point e is a 
discontinuity. The row and column coordinates of e and the continuity count are 
then stored in (IEC(NDIS JC), K = l, 3) . 

4.10 COMraNTS~; 

None 

4.11 LISTINGS 

The listings of this routine, with BOUDIS and JCOUNT are attached at the 

end. 
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4.12 TESTS 


This program has been tested in conjunction with SMOB2, a smoothing ' 
routine documented in the next section. 


no 





ISN 0014 IFn.EU.NREC)NB3«0 









LEVEL 2K8 ( JU)^ 1 


DS/360 FORTRAN H 


DATE 76.159/10.56.29 


ISN..0002.- 
ISN 0003 


111 DJJ AME =L_J1 A IJ^ i Q P1»D 2 » LlNECNtiSS 6 a S I Z E»QQOOKt 

SauRCCfEBCDlC.NDHST.NQOfcCK.LOAD.HAP.NDEOlTt lO.NJXREP 

-SUBROUTINfc .aQUOlSllUaVL.l fciU,Na,NtL»lR.W0li-»lRC*JJ01 

DIHEN5IUN IBDY(NELi3) tIN0nitMB(3} 

kicfmjiariMii.al 


PREVIOUS* CURRENT AND NEXT LINES FOR 'I>1*2|3 RESPECTIVELY. 

-tlJg O-JfH &-D4-SC0N44-KU l - T - .IE S - A T T H E CU R<L€-ig^-4^M&. A-OMC43N4-I WU-I-M 

IS DEFINED AS A BOUNDARY POINT NOT CONNECTED TO AT LEAST TNO OTHER 

DniiAirvAov oniRiTc 


IT IS ASSUMED THAT THE BOUNDARY POINTS IN EACH KOW ARE IN ASCEND- 


!SN 0006 
OShUIiOIU, 
JSN 0008 

I CM nr« 1 


ISN 0011 


ItIpHKiPIKJI 


ISN 0022 


ND2-NB(IND(2H 

JsiaiiNBi-LML-LiU 

IF(N 62 .eQ.O»RETURN 

nn in i=l 


lCDliNT»0 

teller 1 ALJr\ 


» ICQUNT+l 

i t IT R.I D O A \ ' 


, lCOUNT+1 


i;iii jiWBi 





IF( NB3.NE.0) ICOUNT»ICaUNTt 

V . I I Mfi O I -I I «iM«: I » I. I 




ISN 0026 


ISN 0C2B 


Klf i^lllIV^ 


ISN 0034 





IRCINOISilI^IR 




ISN 0030 

1 \W 00?1 1 

IRC(NDiSt3i«ICQUNT 

|0 rriNTINUF - - 


* 

ISN 0032 

RETURN 

M m * M it t /'n 1 IA> Y . ■ t <1 

i it 




wim 




t — ** — ■ P-toi^Y-K-A U * ri :" " " DATS ^7A » 3 ‘8^/14-*r“^^‘SA- 

■ . • . » . 

— tjP'TI Cm'S-—- — ^lAMt-= riAj'^^Tu~P' r ~ 0^ ' f irl*UtrA“NT*^&<>-yS‘f'Z*C'*-04}CrC* K ‘'f ***** ■“ ■ ' 


. SJUKCt iea)»NaUSTi.viuUtCA.LUAD,MrtP.NuEUlT,lD,MaXlitf- . 

iiN .-• j^u.KTi jh jcuoNm x.j»ji-,jz> ooeocec-c — 

ISN i DlMiUSluW 1X(U ,CCOOCOIC 

C" •■ JC ■-■'jfi'I * fiu« 'JF' VrtLUFi {JF ‘JJ SOcH •I HaT—J 1 •L£ ♦ JJ< — — Of 00 CO 30 ~- 

C HtiiJ 1 AtiS( 1 Xt J)-£ Xt JJ) ) . tt. I COClCOAO 

■ - • I i O' OOiT'A — oCiioNl ■(- ■ •■••" - ;•• „ — ; — — •— .i^. •' — —Oi - 

ti'J ! ^F(ol.or.J2>t•gTU^N ' . ’ •* 

ISN K = '‘' ■ — - — - — — — .... — C't 00 C'O 7 O— — 

1)N O'" M DO 10 oHOtOVuC 

ISO O'i-.w - - iruxi jJi~ix( j».Lr-.-i'»oO-ru-K- Oicu-^vc*— 

liH O’! I ' IMl a(JJI-U{J>.OT.UGQ TO 2C '-iCufiCn 

• — •• liN uOJ 3 " -A » 1 — — 1 — C'*'0'j0M C— 

ISN J**!-. 10 CwMINOe OCOCflZC 

. .lii, OOIS 20 JC'-UNT = K ; OvCv/M30 — 

ISN O*'!*!, KLPJffl ’ CC-CC1‘.C 


ISO Of 1 7 ~“tM) — — Z ■ — vCvOOl&O' 




5. SMOOTHING BOUNDARY DATA 

5.1 NAME 
SMOB2 

5.2 PURPOSE 

To patch discontinuities in a digital curve. 

5 . 3 CALLING SEQUENCE 

CALL 8MOB2('lRC, MDIS, IDIS, NDIS, NDEV, lEDY, IWl, IW2, NREC, K) 

where 


IRC is an input array dimensioned (MDIS, 3) with IRC (1, 1), IRC (1, 2)) giving the 
row and column coordinates of the I'th discontinuity and IRC (1, 3) giving its 
connectivity count for 1=1 through NDIS; 

IDIS is the discontinuity number at which the patching should be started (only 
the discontinuities corresponding to I=ID1B through NDIS will be patched); 

NDEV is the logical unit number of a direct access device on which the 
input boundary data set is located; the output after smoothing is written 
back on NDEV. 

IBDY, IWl, IW2 are work arrays to be dimensioned as indicated in the 
listing attached; 

NREC is the number of records in the boundary image; 

K is maximum coordinate difference over which the nearest boundary points 
are checked for patching a discontinuity. (See 5.9, Method). 

All parameters except the work arrays are inputs. 

5.4 INPUT-OUTPUT 

5.4.1 Input 

The input data should be on the direct access unit NDEV, consisting of 
NREC records, the I'th record readable by 

READ(NDEV'I)N, (IBDY(J,1)), J = 1,N). 
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5.4.2 Output 


The output data will be on NDEV in the same format as the input. 

5.4.3 File Storage 
None. 

6.6 EXITS 

No nonstandard exits. 

5.6 USAGE 

The program is in FORTRAN IV and is presently implemented on IBM 360 
using the H compiler. It is available on the user's library in the form of a load 
module. 

5.7 EXTERNAL REFERENCES 

The linkage is indicated in the following table; 


Galling Program 

Programs Called 

SMOB2 

PATCH3 

PATCH3 

SVSCI 

PATCHl 

SORT 

ELIRPT 

PATCHl 

CONTEL 

PATCH4 

PATCH2 

PRTVEC 

SORT 

MVMRMR 

ELIRPT 

VMOV 
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5.8 PERFORMANCE SPECIFICATIONS 


5.8.1 Storage 

The size of SMOB2 is 1068 bytes. Including a main program to supply the 
arrays required to handle a maximum of 2100 boundary points per record with 
K= 20 and the buffers, this program needs approximately 114K bytes for execution. 

5.8.2 Execution Time 

Highly dependent on the image size, complexity and the number of dis- 
continuities to be patched. In the case of the Mobile Bay, Alabama level H GTM 
which h^d 4000 records with 728 discontinuities of which about 530 required patches 
to be generated, the execution time on IBM 360/65 was about 9-minutes. Since 
there is a considerable amount of I/O involved on the direct access unit NDEV, a 
significant improvement in execution time can be achieved by using the array 
read/write routines DARN and DAWN wherever implied DO loops have been used 
in the subroutine PATCH3 . 

5.8.3 I/O Load 
None 

5.8.4 Restrictions 
None 

5.9 METHOD 

The routine SMOB2 simply consists of a DO loop which calls PATCH3 to gen- 
erate the patch points needed for the L'th discontinuity and prints the details of the 
patches produced, with L ranging from IDIS through NDIS, 

Consider the routine PATCH3 . Suppose (I, J) is the address of the discon- 
tinuity at which a patch is to be produced. Then, the records I-K through I+K 
(bounded, of course, by 1 and NREC) are read from NDEV. While each record is 
read one row of a 2K+1 by 2K+1 binary matrix IWl is defined. The elements of 
the row are initially set to 0 and whenever the (J-K+L)'th column in the present 
row of the input image has a boundary point, the (L-rl)st element is set to 1. 

After defining IWl, the routine PATCHl is used to check the array IWl, 
eliminate the I's contiguous with the (K-f-1, K+l)fli element, find fee nearest 1 
among the remaining and join it to that element by a straight line and store the 
row and column coordinates of the points so produced in an array IW2 . Further, 
if the contiguity count of the point of interest is 0, then the I's contiguous with 
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the point joined to the point of interest axe also eHminated and a straight line patch 
is produced to the nearest remaining 1. 

The addresses in IW2 are then merged with the data on the input direct 
access data set by reading the corresponding records of input, sorting the column 
coordinates in each record using SORT, eliminating repetitions of column coor- 
dinates using ELIRPT and writing back on NDEV. 

5.10 COMMENTS 

The routine SMOB2 can be used in conjunction with BOXJDIM or indepen- 
dently. If used independently, the coordinates of discontinuities maybe supplied 
by reading a sequential data set produced by a separate run of BOUDIM. If the 
program terminates due to lack of time, the execution can be continued by a sub- 
sequent run with an updated value of IDIS provided the output data set on NDEV is 
kept. 

5.11 LISTINGS 

Listings of SMOB2 and the important routines called by it are shown at 
the end of this section. 

5.12 TESTS 

This routine has been tested by using the coordinates of the discontinuities 
produced by BOUDIM on the Mobile Bay GTM. The first 40 discontinuities were 
examined in detail by printing the arrays IWl. The performance of ti^ routine 
was found to be satisfactory . 
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LeVFL 21.8 ( JUN I 


05/360 F0RT8AN H 


DATE 76.159/10.57*02 


COHPILER OPTIONS - NAME= HA1N.0PT»02«L1NECNT-56»S IZE»OOOOK, 


.1 

ISN 00P3 
ISW 0004 


SOURCE.tQCDlCtNOLlSTtNOOECK*LOAD»KAP,NOEOir«ID»NOXREF 
SlBROUTl KE 5HOB2(UC.H0tS.lDIS.VtOl 5 .H DEV , IBOY , IM 1 , |H2, NkEC «KI 


COMMON/ PT CHAD/I 1 >J1 .12 «J2.NP 


. iwin».iw2nt 


D'N 16DY(HAX. EXPECTED NO. OF BOUNDARY POINTS IN A LINE AFTER SNQOTHINGl 






LEVEL 

21. B 

1 JUN 74 

) OS/360 fORTRAN H 

"^A te”^. 1 s 57TT. 56 .7T“ 



COMPILER 

OPTIONS - NAME- HA 1 N , 0PT=02 . LI NECMT -56 . S I Z E-OOOOK. 



I .«;w 



SQURCE»EOCDlCtNaLISTrNao£CK,LaADtHAP.NaE01T»lDiNOXAEF 
SlIARnllTINE CONTELI IWl .1M2,H.N>. ‘ 



ISN 

0003 

r 

DIMENSION IHKH.N) «IW2(2tU 





c 

r 

THIS PROGRAM ELIMINATES ALL 1 *S IN IHl CONNECTED TO THE 1 AT 

Ih2 SHn H.E n.IHPNS_inMFn l P,K l UHFRF K ISJUIfF 





C 

r 

THE NUMBER OF NODES IN THE PIECE OF DIGITAL CURVE CONNECTED TO THE 
POINT OF IWTFRFST. .. 


11 

! 

non/^ 

C 

-- - 


.« 

ISN 

I 

0005 

onof. 


L=1 

nn ?n k k=:i ,k . * - 


1} 

ISN 
I 5N 

000 7 
nnoB 


I <>IV|2(1 iKKt 

j5lH2i?.KKI _ . 



ISN 

0009 


iwin»j)^o 


1* 












_L fc VbL 21.8 1 JUN > OS/360 FURTRAM H OAT£ 76. 15 9/lQ« 56. 40 



CGMPILER 

OPTIONS - NAME* HA! N •0PT»02 . L INECNT »56t S I Z E-OOO.OK, 



liM_ 

non;^ 

SOURCEf EBCDlCtMQLI5TfNOO£CK«tOADiHAPfNOEDITf I0« NQXREF 
SnflRnUT ! MF PATfRl < f T U-? .fM tWf I Tf nUMT 1 



I SN 

0003 

COMMON/PTCHAD/Il tJlf I2,J2»NP 




J1C0A 

OAXA-J-EJlS^ZOi! 




ISN 0005 1PASS«1PASS+1 

■l.S li - 0 Q04) l-aUP4.S-4^i. E . A QKAL L P ft TVE C ( W -tWm-,44 

C 

C (2<‘1CUUNTI NDNCONTICUOUS NEIGHBORS. 

C 

C SEE CONTEL FOR DIMENSIONING INPO FOR 1U2. 

C 

ISN 0008 DlHENStON I HI (M .Nl » I W2 ( 2 .1 > 

4-SN-OOOa l4t2-(4-4-»-a4 

ISN 0010 IH2(2»l)3j 

C 

C ELIMINATE POINTS CONTIGUOUS WITH U,J). 

C 

ISN 0011 CALL CONTELUUl ,1W2»H,NI 



ISN 0023 1F( Il.EQ. OIRETURN 


C PRODUCE PATCH ADDRESSES IN IH2. 








LEVEl. 21. B ( JUN 74 I OS/360 FORTRAN H DATE 76« 1 S9/10» S6» 


COHPlLEft OPTIDWS - NAHE» MAI t4,OPT «02 . H MECNT »56 , S 1 2£«OOOOK, 

5CUPCEiEBCOICtNQLISTfNOOECK>LUAO»HAP»NaEOIT»lD»NJXKEE 

I SN 0 002 SiJaSMLLN£-.P-AlCtl2iJ-tLtllP.>.l 1 1 1 1 2,J2i : 

iSN 0003 DIMENSION IW(2.U ! 




c • 

r: 

TO GENERATE COORDINATES OF LINE JOINING IfltJll AND (12>J2). 
niMFNSION I WI_2^,MP1 HHERE HP* MAX._NO._ OF FOINTS EXPECTED TO BE 



c 

c 

PRODUCED. NP- NO. OF POINTS ACTUALLY PRQUUCEOiBV THIS ROUTINE. 

ISN 

000^ 
0C0 5 


IHLsMINOni tI2» *1 

IMX=MAXO( 11 .12»-1 

r SN 

0006 


JMN^HINOI J1 « J2K1 



ISN 0026 IH(1»NPI>1 










LEVEL 2U8 ( JUN 7A I 


DS/360 FORTJIAN H 


DATE 76. 1S9/10. 56.49 


COMPILER OPTIONS - HAHE= HAl N. 0PT»02 » L INE CNT 1 S I Z E» Q Q0QKt . 


. LSM...OCfl2_ 


!SN 0003 


SOURCE tE8CDlC,N0LtST»MQ0ECKflUAUi MAPf MOLD IT« { Dt NOXR EF 
_SiJllRUUl.I.NE._E.AJa^3-(j-aOl^,J.^J.,J.CiJUi^L,^C^i^*Jl^^il2,JJR£C,J^0£«J 


COMMON/ PTCH AD/n, Jl ,12. J2.NP 
.n.i.ME 


OF P ATCH P-OUiLS- 


EXPECTEO TO BE GENERATED BY PATCH2, DIMENSION REOUIREO BY CUNTEU 
PATHH fltSrnMTIMllMV 


ISN 

-LSN. 


0004 

-DOOS- 


ISN 
r <iM 


0C06 
nr, 07 


DIMENSION iwllt I.Iri2(2.1l,IBDY(U 


Kl=HAXOn-K»lI 
K 2 = M INO ( Lj- IL. NBh CJ^ 



ISN 0021 IMimBDYlLI-lCLHMKj*KK2UKK-Kl*n»l 

LSN-00Z2 ZO CaNXIiUlE 

K ISN 0023 10 CONTINUE 

C 

C GENERATE PATCH ADDRESSES IN IH2. 

C — 

ISN D024 CALL PATC HU INI .1 N2 ,KK2 1 ,K2 1 ,I-Kl ♦! ,K»l , KOUNT I 
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LEVEL 21.8 I JUN 74 » 


U5/360 FOKTRAN H 


DATE 76.159/10.56.5S 



ISN 0013 


ID»( 











HBHIi 

1 

1 

1 

1 

1 

1 

1 

"TaTE^T,1597!T56,?^ 


L£Vht 

21. B ( 

JUN 74 

t U5/360 FORTRAN H 


COMPILER QPTIOMS - NAME- MA I N , OPT «02 , L I NE CNT «5 6, S 1 Z E'OOOOKj ! ^ 

1 SN 

000 2 

. 

SOURCE.EBCOIC ,NQLI ST t NJOcCRt LU ADt HAP,NUEO ITi I O.NJXREF . 
•JUnRnUTINR PRTVECI IX.N. IFHT > - 



ISN 
I 5M 

0003 

0004 


DIMENSION IXINI 

IF( IFMT.FQ. 1 1 PRINT 1 00 . I X 



I5N 
I .SM 

0006 
000 8 

• 

IFlIFHT.EU.aiPKlNT 2C0,IX 

RFTimW . . 



ISN 
I SM 

GC09 

0010 

100 

200 

FQRMATI 10X41111 
FORMAT! 1X4013) 



ISN 

non 


£N0 ^ 





6. IDENTIFICATION OF CONNECTED REGIONS 

6.1 NAME 
REGIONS 

6.2 PURPOSE 

To identify all distinct connected regions in an image given the boundary 
data in SLIC format and produce a map \vith a number at each point showing the 
region to which it belongs. The region numbers will be in descending order of 
area. 

6 . 3 CALLING SEQUENCE 

This is a main program. In its present version the image size is supplied 
through DATA statements. 

6.4 INPUT-OUTPUT 

6.4.1 Input 

The input to this program is a sequential data set on lo^cal unit 8, having 
NREC records stored as N, (IX(J), J = 1,N) in unformatted FORTRAN mode. 

6.4.2 Output 

The output of this program will be a sequential data set on logical unit 12, 
having NREC records with NEL pixels each, with one half-word (2 bytes) per pixel . 

6.4.3 File Storage 

This program requires a direct access data set with NREC records and 
NEL half-words per record. 

6.5 EXITS 

Not applicable 

6.6 USAGE 

This program is in FORTRAN IV and is implemented on IBM 360 with the H 
compiler. The associated subroutines are available as load modules on the user's 
library . The deck for the main program is available with the authors and needs 
only slight modifications in the DEFINE FILE and DATA statements for use on 
any data set. 
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6 . 7 EXTERNAL INTERFACES 


This program uses several subroutines as indicated by the linkage table 
below: 


Calling Program 

Programs Called 

REGIONS 

PET 

VMAXI4 

VMINI4 

RIDER 

SVSCI 

DARN 

SEQLS 

SAWN 

RIDER 

SVSCI2 

SVSCLl 

SORT 

RIDERl 

RIDER4 

DAWN 

VMOV2 

RIDER2 

PRTVE2 

DARN 

VMAXI2 

SEQLS 

SORT 

FIIPV 

SORT 

MVIVIRMR 

RIDERl 

SVSCI2 

RIDER4 

RIDERS 

SVSCI2 

PRTVE2 

RIDER7 

— . 

RIDER6 

SVSCLl 
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Calling Program 

Programs Called 

RIDER2 

PRTVE2 

RIDER7 

VMAXI2 

VMrNI2 


6.8 PERFORMANCE SPECIFICATIONS 

6.8.1 Storage 

The present version of the main program is 134, 436 bytes long. The ex- 
ternal references required and the buffers increase this to 192K bytes . However, 
the- size is dependent on the data set to be handled and the dimension statements 
should be changed to satisfy specific requirements . 

’ DIMENSION , 1^" *lip;s‘(MSEG+l).J~’ 

INTEGER*2 IW1(NEL), ’lW2(NEL), iTABL(MR*MSEG), IS(MR) 
INTEGER*2 LW(MR) 

LOGICAL*! rDENT,(MR,MR) 

where 

NR ' Maximi^ number-^ region expecte d; 

N’_ = JVIaximum number of boundary points expected in a record; 

MSEG = Maximum number of "segments" required to handle the 
image (see 6.9, Method); 

MR = Maximum number of regions identifiers permissible in a 
segment; 

NEL = Number of pixels per line in the output map. 


6.8.2 Execution Time 

The time is highly dependent on the size and complexity of the image. The 
Mobile Bay GTM (level II) resulting in a region identification map with 400x2100 
pixels and consisting of 17 42 regions had to be handled in 15 sections and took 
19.5 minutes of CPU time on IBM 360/65. 

6.8.3 Restrictions 
None 
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6.9 


METHOD 


This program has five major sections. 

(i) Determinatton of the bounds on the column coordinates of 
boundaries on the input data set; 

(ii) Finding a preliminary set of region identifiers; 

(iii) Finding the areas of each of the regions; 

(iv) Generating a mapping such that the region numbers are used 
in the order of decreasing areas; 

(V) Modifying the region numbers by table look-up. 

6.9.1 Determination of Bounds 

The maximum and minimum values of the column coordinates of the boun- 
dary points are determined. If the minimum is greater t han l, it is set to 1. If 
the. maximum is less than the value of NEL supplied, it is set to NEL. The value 
of NEL is then changed to Max-Min+1, The output image size will then be NREC 
by NEL. 

6.9.2 Finding Preliminary Region Identifiers 

This is the most important step in the program. The subroutine RIDER 
is used for this purpose. Its function is similar to -the routine with the same 
name described in [4] , The routine in [4] was designed to print an error 
message and return with NR=0 when the number of distinct regions exceeded 
MR. But the present version can handle up to MR*MSEG distinct regions while 
still using a "region identity ’matrix" of size MR by MR (rather than MR*MSEG 
by MR*MSEG). 

This routine uses the arrays IWl and IW2 as the previous and current 
records of region identifiers . By convention, region numbers 1 and 0 indicate 
the "exterior" of the image and boundary points. The MR by MR array IDENT 
is used to store information about identify of regions, IDENT (I, J) = .TRUE, 
meaning that region numbers I and J refer to the. same connected region. 

Initially, the array IWi is set to all I's and IDENT is set to all .FALSE. . 
Each of the input records is read and the following operations are performed. 

The boundary coordinates in the input record are arranged in ascending 
order. The routine RIDERl is used to generate, in IW2, the region identification 
numbers corresponding to the present row. First, all the elements of IW2 
corresponding to the boundary coordtnates are set to zero; Each interval between 
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the zeros is compared -with the corresponding segment of IWl. If there is no non- 
zero element in that segment of IWl, a new region number is started and assigned 
to the interval in IW2. If there is a nonzero element, that number is filled into 
all elements in the interval. Finally, IDENT(IW1(I), IW2(I)) is set to .TRUE, for 
1 = 1, NEL wherever IW1(I) 4 0 and IW2(I) 7 ^ 0, indicating that IWl(I) and IW 2 (I) 
refer to the same region. Also, when new region identifiers are to be used, the 
routine RIDERl verifies whether the number of identifiers exceeds MR. If so, 
the value of NR, the total number identifiers, is set to-NRP, the total number 
up to the previous record and the control goes back to the routine RIDER. 

Now, if RIDERl returns a positive NR, the array IW2 is written as the 
I*th record on the direct access data set (unit number IDEVO in RIDER, same as 
90 in the main program) and IW2 is moved into IWl (so that it becomes the 
"previous” record while handling the next record) . 

If RIDERl returns a negative NR, then NR is changed to -NR and the routine 
RIDER4 is called. The set of records handled between any two calls of EIDER4 
will be referred to as a segment. Associated with each segment, a table is defined 
which gives a mapping from the set of region identifiers obtained in that segment 
to a new set reflecting the connectivities discovered up to the most recent segment 
handled. Also, the initial record number for each of the segments is stored in an 
array . The functions of the routine RIDER4 are to: 

(i) Reduce the matrix IDENT (using RIDERS) examining all of the 
available connectivity information in it and obtain a look-up table 
for the current segment; 

(ii) ’ Modi;^ the t ables for^thie previous segments to reflect the newly 

jfojund cormectivitiesv, if_any; - 

(iii) Find all the distinct -re^on numbers occuring in the last record - 
IWl of the current segment and change the numbers there which 
are greater than 1 to consecutive numbers starting with 2; Let 
NR be the largest number in IWl; 

(iv) Set up an array IS consisting of the distinct region numbers in 
IWl and then change IS(I) to ITABL(IS(I), ISEG) where ITABL 
is the look-up table for .the current segment; 

(v) Set ail elements of IDENT TO .FALSE, except when IS(I)=IS(J) 
for I, J in the range 1 through NR. 

After each call to RIDER4, the segment count ISEG is incremented and 
the initial record number for the next segment (which is really the record number 
at which RIDER4 had to be called) is stored in IRES(ISEG). If MSEG is exceeded 
by ISEG or if NR> MR (which means there are more than MR distinct regions in 
the last record) the routine RIDER prints an error message, sets NR=0 and exits. 
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otherwise, RIDERl is called again, IW2 is found and written on IDEVO’ and the 
pragram p:raceed^ no^rm^y to_&e nextjnput_re^rd. __ 

After the NREC input records have been processed the routine RIDER4 is 
called to get the look-up table for the final segment. A call to RIDER2 changes 
the look-up tables for all the segments such that consecutive region numbers are 
used. 


Finally, each record from IDEVO is read, the appropriate look-up table 
is used to modify it and the record is written back on IDEVO. Also, NR is set 
to the maximum region number used after table look-up. 

6.9.3 Finding Areas 

A histogram of the region identification maps is found, giving the total 
number of occurrences of each of the region identifiers 0 through NR. These 
numbers indicate the areas of the regions. 

6. 9.4 Finding the Final Look-up Table 

A sequence of natural numbers is used as a secondary array with the. histo- 
gram as the primary array in a descending sort operation (routine SEQLS) . The 
resulting secondary array then gives the sequence of original region identifiers 
corresponding to decreasing areas . An inverse mapping [inverse mapping of 
{K(J) J = l, N] is defined as (IY(J) J = 1,N] if lY (IX(J))=J.] of this sequence gives 
^e final look-up table. The actual coding follows these principles' but is slightly 
difierent in detail to preserve the identities of regions 0 and 1 which have special 
significance. 

6.9.5 Deriving the Final Re^on Identification Map 

The look-up table generated above is used to modify the region identifiers 
on IDEVO, record by record', and write out the final sequential data set on unit 12. 

6.10 COMMENTS 

An approach suggested in [5] can be used instead of the one described above. 
With that method, the processing would be identical, except that the matrix IDENT 
is not defined. Instead, a table is updated every time a new connectivity is dis- 
covered. While this saves storage, it appears to take more execution time than 
the present method. 

6. 11 LISTINGS 

The listings of the main program and the associated routines are attached 
at the end of this section. 


131 



6.12 TESTS 


This program has been tested on the Mobile Bay GTM both, before and after 
smoothing and found to work satisfactorily. Also, the results have been found to 
be identical (on a smaller data set) with those obtained by the earlier version of 
this program. 
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as/360 FURTRAN H 


LGVtL 21*7 ( JAN 73 \ 


UATE 76. 163/U. 53*57 


COMPILER OPTIONS - NAME= HA I N tUPT*00 t LINECNT »56 r 5 Ut»OOOOK# 


GO 

CO 


ISN 
I SN 
ISN 
ISN 
ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
I 5N 
15N 
ISN 

ISN 
ISN 
I SN 
ISN 
! SN 
ISN 
ISN 
ISN 
ISN 
I SN 
ISN 
ISN 
ISN 
ISN 
I SN 


ocoa 

CC03 

0C04 

nC05 
CCD6,„ 


SGURCEtEBCDIC tLIST»N00£CiVfLGA[)rHAPtNUEDn«lDiN3XREF 

.DIMENSION I XI4000I .IRESian .. . . 

INTEGER ^2 I WU 2 1 00 ) 1 1 ^2 ( 2100 I fl TA6LI 8300 I » IS C400I 

.INTEGER *2 LWUOO) 

.LOGlCAL^l IDENT(300,300) 

-D££INfL..FJ LE .9CI <^0U0i^20Qf L# I AVI 


.aiN. IXlMAX(2NR^2fN) HHEKE NR^HAX* NO. OF REGIUNS^EXP-ECT £0- AND . 
N* MAX NO. OF BOUNDARY POINTS EXPECTED IN ANY KECQROI 


ono7 

ocrs. 

C009 

COIO. 

coil 

CC12 

cm 3 

0C14. 

0015 

CC17. 

0018 

C019, . 

0020 

oczr 

C022 

0023 

002^ 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 , 

0033 

0034 

0035 


DATA NREC •MRfMSEG/4000,300,20/ 
J)ATA.NEL/2100/ 


aiL 


6GQ„ 


HAXX=-10C0000 

...minxt 1 cococo 

CALL 5TRTHR 
..CALL PETIO) 

DU 10 I =1 ,NR£C 

.READieiN, nxui ,J«i,N) 

IFIN.EU.OICO TO 10 

.CALL VMAX14nXiN»HAXXl 

CALL VM1NI4( I Af NfHINX) 

..CONTINUE , 

REFUND 8 

.IDENTIFY CONNECTED REGIONS^. 
PRINT 600fMINX,MAXXfNEL 
.M1NX=MIN0 (MINXp1) 

PRINT 600 ,HI NX,MAXX,NEL 
.MAXX=MAXO(MAXXrNEU 
PRINT 600tM!NXrMAXXtNEL 

.NEL^HAXX-KINX^l 

PRINT 600 tMINXf MAXXiNEL 


ISN 0036 


...FQRMATC* MINX,MAXX#N£L»^3ia I . . 

PRINT lOOfNRECfNEL 

..aCJ} FORMAT!//* IMAGE S UE =« I • 1 5 , • • ! I S • M • I 

N0UM=1 

PRINT lOOOiNDUrt 

FORMAT! • NOUM-* 151 
.CALL PET(2> 

CALL R1 DER( 8 f NR EC «NELt9Ct Ml NX «IX» INI tl»2 I IT A6L« IDEN 

. MSEC rlRESfISI . . 

CALL PETI2) 


1000 




M 


T I MR ff NR f LN t 




FIND AND PRINT HISTOORAh OF REGION IDENTIFICATION HAP< 


dr 


1 SN 
ISN 
1 SN 
ISN 
ISN 
I SN 
ISN 
ISN 
ISN 


0037 

0038 

0039 
C040. 
C041 
C042 
00^3 

0044 

0045 


PRINT 2C0 

2C0 F0RMAT!//10X*REGI0N NO. * lOX * NO.^OF PIXEL5M*-- 

CALL SVSC I I IXfNR^l ,01 

D0..30..I =1 ,NREC 

CALL DARN!90,1,IW1 ,NEL<»2I 

DQ.30 IEL«1,NEL . ... . . 

J^IWU lELUl 

IX(.Ji»IX! J>+1 ... ... ... 

NR1=*NR^1 




1-d ^ 


.J_Q 


PAUe 002 


!SN C0A6_ 
ISN P047 
ISN C040 
ISN 0050 
ISN 0051 
ISN 0052 


AO 

-.300- 


_Dli. A0,I*1,NR1 

J»I-1 

_.1P( ixn t. NE. 01 PRINT 300.Jilxm 
CONTINUE 

- FORMATinXIA, 16X191 
CALL PET(2I 


REARRANGE NUMBERS IN DESCENDING ORDER OP POPULATIONS. 

-LEAVE 0 AND 1 UNCHANGED SINCE THEY CORRESPONO-TQ-EXTERIOR AND — 
BOUNDARY POINTS RESPECTIVELY. 


ISN 0053 CALL SE ULS( I X( 3 1 f I XI NRl^l ) »NR>I iNK-II 

„ISN ODSA PRINT. AOO 

ISN 0055 AOO PORMATI'l REGIONS AFTER REASS IGNHENTS: • I 

, ISN 0056 PRINT 200 • 

ISN 0057 00 50 I tNRl 

ISN 0058 JsI-1 — ... 

ISN 0rS9 IFI I.LE.2IPRINT 3C0*JfIXUI 

-ISN 0061 IFII.GT.-2 >PRINT-350^J^X< U-r4»44>NR- l > 

ISN 0063 50 CONTINUE 

ISN 006A 350 FOSMATl ll XI 6.16X19.17X161 

ISN 0065 ' 1X11)1=0 

ISN 0066 IX(2)-l . . 

ISN 0067 DO 60 I -3.NR1 

-.1 SN .0068 60 ULLXIitNR-:!.! <-31 

ISN 0069 CALL PET( 21 


ISN 0070 
ISN 0071 
ISN 0072 
ISN 0073 
ISN 0074 
ISN 0075 
ISN 0076 
ISN 0C77. 
ISN 0073 


80 

... 70 — 


MODIFY REGION NUMBERS ACCORDING TO NEN ASSIGNMENTS FOUND IN IX. 
00 70 I =1 ,NREC 

CALL DARN(90.I»IU1,NEL«2) • 

DO 80 lELsl.NEL 

IUinELt=>IX(J) 

.CALL-SAHNa2.IWl-tNEL*2» — I- ■ 

CALL PETI2) 

END 



?si: 











t fcv cL /-i—jnti -79--J tiS/3«9 — FtiKTft A^--H ; OA-fd 

... MwMir* — M£CNF-ifrrS-J 2t-»0&eC-H-f ; ! — 

St)LIKCtfoCU,NULI5r,UuDECK.LUADfHAl'*NCit;tJlTtIU»>jaXk6f l . . 

"*■* ■* I Srt u 111 Iiijfc — k V^Ik ( a t H A*y-lH f & y^t b»lA»IB-^ — — ■ ■■■ . — i ■ — *.». — ■ — — — — ■ — . — _ 

isij :'(C>3 uri'cfoiou A(MA,m tUnutNi .... . 

ISM Ov-'A Oil Ij.- 

1 00' S Odo , J>-A( lA.JI 

....... 15f4 |.tIor<W ; 

ISO CC{ 7 • Enii 




SUBROUTINE PETd) i 

j F <- 1 .. £ .-0-) GU- TO -1‘0 

CALL TlHERdTlHEU 

TTIHE = 0. 

WKirE(6,20CJ 

RETURN 

■ — 10 — CALL rfH£RUT!KE2) 

TlMb*( IT1HE2-ITIHEU/100, 

TIlrtt^IIlME+TIME 

ITIME1-IT5HE2 

KRi rt('6— it'cmtiE-rt-rt-ME 

100 FURfUK lOX'TlHE ELAPiEO SINCE LAST PRINTINO Of TIHe»*£l2.3. 

— lOTAL-TlME -eL-APStO.-<El-2T3-,-*5-eC— 

RETURN 





> o 



e — 

c 

f* 

1 L MN ♦ L c « Mi fn-i-TJ U ¥ -ru J T i ♦ L t: ^ X t « Ni: l— i" 

DEFINE FILE 1 DE VO ( NREC ,NE L*2 , L, I A V I 

C 

DIMENSION IBOYII) 

1 nr T r- A 1 r_n jrMjr u Q \... 

- 

tOGICAL 1 LW(.MR,MR) 

T t ).. T W7 { mp-i-wt TABt r 


DIMENSION IRESIMSEG) 


INITIALIZE A WORK ARRAY IWi WITH I’S AND IDENT WITH .FALSE. 


CALL SVSCI2nWl,NELrU 



LOOP ON RECORDS 


DO 10 IREC=1,NR£C 


C READ ONE RECORD OF BOUNDARY INFO. 

■e : 

REAO(NTAPB)Nr {I8DY( I ) ,1=1 ,N) 

If - I - NvNEyQ TC ALL SORT ( I DDY, 1 ,N , N ,1 ,T ,TT ) ^ 

C 

€ U - S€-TW1 AND IDDY TO SET ARRAY IW2 AND MATRIX IDENT 

% 

C 

- 3 ^ ^ 

CALL RIDERl ( IWltlBOY, ICMN,NEL ,N,IW2,I0ENT,MR,NR}- 

I r (NR .GT. - 0)"G0"T0 2 0 

PRINT 2C0, IRECtNR 

F -e RMAT (» (IREC - yNR) = <216 1 

IFIIREC.EQ. IRES! ISEGHGO TO 40 

^ 

CALL RIDSR4{IDENT,LW,MR,NR,ITABL, I5EG, IWI ,NEL , IS , .FALSE . ) 

I -- 5EG=I5EG -<‘ 1 

IF( I5EG.GT. MSEGIGO TO 

ITl-K- H SF6-H^T-R eC ^ 












IF(NR,LE, MR )G0 TQ 30 


Hv: r irr r i i — rt i r> l k; 

NR = 0 ) 

ir&O p-0RWH-+—tRic8R— Ca Na-I“Ti-&N— i-N— R-t D-ERi -SUi*?" !rI-E-&-MR— eX'CE c £ 

.EC AT RECORD NUM 86 R‘I 6 ,«. RETURNING WITH NR= 0*1 

R-ETOR-fJ 

20 CALL DAWN( IDEV0,IREC,IW2,NEL*2> 

e - A t- L- VH -& V2( ’ ir2 ' ,NEL,Iwl ) ■ 

10 CONTINUE ' 

e-Mr lr H ^T -^iRV<-t-&£iiTtL W f MR tN Rf KA r6tri-5£-G^ IHl tNELt ISt . - TR -UE^ 

CALL RIDER 2 (lTA 8 L»MR<'ISeG) 

^frE^- ( T 51 -G * 1 ) ^ ' 

PRINT 3 G 0 

-3^-5 FORMA-T(/ - // - ' - F - I - f )- AL - - TADLE - ^ '- fOR MSDI - FYIUG Re&TGN - NU - MG r R - S -*-) 

DO 60 jSEGsltlSEG 

P R i N - T - ^C - 0 -- vI - 5 E -6 — — 

AOO FORMAT! ‘OSEGHENT NUMBER* 13) 

-irO ^Atir-PR^-E-2i-i-T^ { 1 1 J 5 E G — 

JSEG=1 


UX) ^ -XL yT'f uX 

IF( IREC.LT. IRES! JSEG+inGO TQ 50 






FC^TIUI. II 


tEVCt Ti;7' JAN- 7? ) 

/ 

eirir JiLr:’t'Tji^rtin'rl^^TJAnir= — nn iirt l i r^rtm-r 6“#5 T iM-vr: : k r 


"OTVTE”'’ 7 crn‘ 27 tnr;T 97 Tr 


— 'C ^2 


- ,JSN 
ISN 
IbU 
ISN 
iSN 
l:»N 
“iSU 
15 N 


(h)c 3 -- 

0 CO 5 ' 

T'UOV 

cei:*' 

00 n 


S u Ur^C t « bC J t ulI L I S T f Muut C K f L c i^D t H •At' t •% LtO 1 T 1 1 0 1 Wux K £F ’ 
— iu:?KutJ 1 UJc'-KtoiTKi t iiCf IdOY , IttiNvNtL# kf I Y'f lutrM irtKfHK) 


1. 1 Vt fr tUHKCi\r' it ? ' tJt* bJUNOAHY ADUK’tSiti ( i U UY U ) # I - 1 f N) * ANU THl; 

lA:t\ A^^^AY IXt r i CUK^tkf AivkAY IY CuNTmIMAw rtejlUNi IUlNTIFICA-' 
Tlj»j wurtncKi* *A4.'iur IP f«c kukeuirkUAKt tLt«tnli liv cukkcnT' kuw — 
akI; CljiiriuLlJUS <«ITh At«Y Uuito LOuDAKV PulNli ut Tnt; LAiT rtOhf 5 ^T 

CuK^t it'l'kUl Nu tLcHcwl^ IN tptr«i hnWiX. ‘ - -- 


* 1 - 


:oi2 

OCl 0 - 


c 
c 

“C — 

c 

C"' Int 

C 

bir-,t Ni lutt^lou Y(N> « ^ , ' ■ 

luTLbi:^ IXiNwL) pi Yl ;y£U 

- CU^ilCAL^l lafcKUUKpnK) ‘ ““ * 

* lriN*k(:ic)ob Tu 10 

r p"P/t L tr> 

Cj fu 2 ‘ 

IC^ — C-jr.T IwUt — ^ — — — 

Urv P - NK 

C ALL PuInT5 Fu Frtfc LtP7 Of JflDYll) AR£ •tXTcKlLK^ Pjl/vTS iHtOWU 1) 
1 = IduY( 1)-ICHN 

iM i.oTW’ KALu-^v*5Ci2nYf n rr 


t 5 **r 

Kisn 

.. to — - 


iSN 

liN 

15 N 


3016 


C018 

O^'lV 

C020 


"C Axu'poirsry ra“TK£'^'R:rGHr‘ar t^oy nrr *' ake' »EXTE;ncjK» ‘rciHrs'r 

c 

|^luDY(«)-ICMN*2^^ * ' 

IPlNcLsod. nCALL SViCI2UYm pNLL-l<-i fl \ 


" C“" 

c 

-c- 


3C 

' C - 

c 

. ^ . 
c 

'*C*‘ ’ 


DtSiOMArc ALL JUUNOA^^Y FUiiillS AS *K£:bIDN P*. 
DJ 30 l=!pN ■ 

'J=I»0YU J- ICMN^l 

lYU l=^C 


HjR I = lpN-l fcXAMlkt IX(J) FUR IdDYC IKlT.J^LK XdUY( I + U 

mO‘StJ lY ACLUKDINLLY* TJfW“KFunrN TYtJtiE LK' 1 5 'S f AKTtO ’ WH£N“ 
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CALL PRTVE2(ITABLC1,2>»NR> 


CHANGE ITABL(*,1):. 


DO 30 I=1»NR 


RETURN 
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ING (IREC + IPTH RECORD EXCEEDS MR. { L AST =. F A LS E . ) . 
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NCR=NR 


CALL PRTVE2{IS,MR) 


I S( I )=ITA8L( IS( U f ISEG) 


CALL PRTVE2USfNRI 


SECTION 4 


BY TESTING WHETHER I S U ) . EQ. I S ( J ) 
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L U ft. .* X X f IN 
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iW 
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W VJ A A, f IN*N 
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RETURN 
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7. DELETION OF BOUNDARY POINTS 

7.1 NAME 
DBOUND 

7.2 PURPOSE 

To modify each, of the ’'O'* pixels in an image to the most frequently occurring 
number in its 3 by 3 nei^borhood, (This is useful, for example, in generating a 
level I GTM from a level II map and/or suppressing all the boundary points in a 
GTM and replacing them with reasonable class labels). 

7.3 CALLING SEQUENCE 

CALL DBOUND (NREC, NEL, NEL2, NTAPI, NTAPO, IK, lY) 

where 


NREC = Number of records in the input image; 

NEL = Number of pixels per record; 

NEL2 =NEL+2; 

NTAPI, NTAPO are the logical unit numbers of input and output sequential 
data sets; 

IX, lY are work arrays to be dimensioned as indicated in the listings. 

All the calling arguments except IX and lY are inputs. 

7.4 INPUT-OUTPUT 

Both the input and output sequential data sets have the same format. The 
number of records is NREC . The number of pixels per record is NEL and the 
number of bytes per pixel is 4. The records are in unformatted FORTRAN. 

7.5 ■ EXITS 

No nonstandard exits 

7.6 USAGE 

The program is in FORTRAN IV and implemented on the IBM 360 using the 
H compiler. The program is in the users' library as a load module. 
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7 . 7 EXTERNAL INTERFACES 

The subprograms required by this routine are: 

SARN, a sequential access array read routine; 

VMOV, a routine to move a vector in core; 

MAJOR, a function giving the most frequently oc curing number in a 3 by 3 
nei^borhood. 

7.8 PERFORMANCE SPECIFICATIONS 

7.8.1 Storage 

This subroutine is 1036 b 3 d;es long. With the main program needed to 
call it for an image with NEL = 866, the external references and buffers, the 
storage required is 40K bytes . 

7.8.2 Execution Time 

Depends on image size. For a test case of 1624 by 866 pixels it took 
approximately 100 seconds. 

7.8.3 I/O Load 
None 

7.-8. 4 Restrictions 
None 

7.9 METHOD 

This program uses a circular buffer IX with pointers II, 12, 13 indicating 
the previous, present and next records under consideration. Imtially, II, 12, 13 
are set at 1, 2 and 3 respectively. After each record is processed, the pointers 
II, 12, 13 are '’rolled” upward. The processing of each record consists of checking 
the eight neighbors of each pixel whose value is zero . The function subprogram 
MAJOR is employed to determine the most frequent number occuring in the set 
of eight (If such a number is not unique, the first encountered number is taken) . 

Records -0 and NREC+1 are defined to be identical to records 1 and‘NREG""®=®^ 
respectively. Also, pixels 0 and NEL+1 in any record are defined to be the same 
as pixels 1 and NEL in the same record. 
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7.10 


COMMENTS 


None 

7.11 LISTINGS 

The listings of DBOUND and MAJOR are attached at the end of this section. 

7.12 TESTS 

This program was used in removing the extraneous boundary points after 
conversion of the level n GTM of the Mobile Bay region to a level I map. 
Line-printer plots of the maps before and after the application of DBOUND in- 
dicate satisfactory operation of this program. 
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8, COMPUTATION OF CONTINGENCY MATRICES 

8.1 NAME 
CONTMATS 

8.2 PURPOSE 

To obtaia and print "contingency matrices'! 6] , showing, for all pairs of 
classes in two classification maps, the numbers of simultaneous occurences of 
various types of transitions (no boundary, horizontal boundary, vertical boun- 
dary and boundaries in both directions) . 

8.3 CALLING SEQUENCE 

This is a main program which takes card inputs as indicated below: 

READ 100, TITLE 
READ 200, NREC, NEL, M, N 
100 FORMAT (80A1) 

200 FORMAT (416) 

where 


TITLE is a title of up to 80 characters to be printed on the top of 
each page of output. 

NREC, NEL are number, of records and number of pixels per record, 
respectively. 

M, N are the numbers of classes in the two maps. 

8.4 INPUT-OUTPUT 
8.4.1 Input 

The input maps should be sequential data sets on units 8 and 9 with NREC 
records and NEL pixels per record on each of them. The number of bytes per 
pixel should be 4 and the records should be unformatted and FORTRAN readable . 

8.4.3 File Storage 

None 


B.' 

01 
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8.5 


EXITS 


Not Applicable 

8.6 USAGE 

The program is in EORTRAN IV and implemente.(Lon the IBM 360 using the 
H compiler. The program and the external references needed are in the users'' 
library as load modules. 

8 . 7 EXTERNAL INTERFACES 

The subroutine linkage is indicated in the following table: 


Calling Program- 

Programs Called 

CONTMATS 

CONMAT 


CONMXP 

CONMAT 

SARN 


VMOV 


SVSCI 


INTLOG 

CONMXP 

SVSCI 


8 . 8 PERFORMANCE SPECIFICATIONS 

8.8.1 Storage 

The present version of the program is designed for NEL ^ 2000, 15, 

Ns 15. The main program size is 62340 bytes. Including the buffers and the ex- 
ternal references, the program requires IlOK bytes. 

8.8.2 Execution Time 

Depends largely on image size and number of classes . For the case 
NREC = 1624, NEL = 866, M=N=6, it took approximately 10 minutes. 

8.8.3 I/O Load 
None 
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8.8.4 Restrictions 


None 

8.9’ METHOD 

The definitions of contingency matrices used here have been discussed in 
detail in [63 and will not be covered here. The subroutine CONMAT is used to 
find a four dimensional array EH [di mens i oned ( 4, 8, M, N>] and the routine 
CONMXP is used to print 

(i) M*N matrices (size 4 by 8) showing counts of agreements and 
disagreements for each type of transition for each pair of classes; 

(ii) M*N matrices (size 4 by '4) showing counts of each type of 
transition for each pair of classes obtained by adding the 
right and left halves of the corresponding matrices from 
(i) and dividing by .3; 

(iii) A 4 by 4 matrix showing totals of each type of transition 
obtained by adding all the matrices in (ii) ; 

(iv) An M by N matrix which is the joint histogram (contingency 
table) of the two input maps, whose (I, J)th element is obtained 
by adding all the 16 elements in the 4x4 matrix corresponding 
to classes (I, J) defined in (ii); 

(v) The individual histograms (inventories) of the two maps ob- 
tained by adding the columns (for map 1) and rows (for map 2) . 

Also, the transition and point by point similamty counts (traces of the 4 by 
4 and M by N matrices, respectively) and percentage similarity measures are 
printed. 

8.10 COMMENTS ' . 

None 

8.11 LISTINGS 

The listings of CONTMATS, CONMAT, CONMXP and INTLOG are attached 
at the end of this section. 

8.12 TESTS 

Portions of the printout from a test run on the Mobile Bay GTM v/s LCM 
(December 5, 1973) are shown at the end of this section. 
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REAB( 5, ICC »TITL6 





C 

c 

— — ^ — ' — * 1 ‘ ' 

FIND AND PRINT MATRICES SHOkING NbPBFRS UP AGR.EEMNETS AND DIS- 




15N 

.CJ06 

c 

AGREEMENTS OF TRANSITIONS FOR EACH PAIR OF CLAjSSES. 
CALL CONH4T(3,9,NREC,NEL.M.N. IX.IY.IH) 



w 

ISN 

I6N 

0)07 

0038 


CALL CGNHXP(riTLE,IH,IH2,lNVl,UV2,NREC ,NEL,M,N) 
STOP 1 



V) 

M 





LCVfci. ai.9 


isKjjoa 
‘isN dooi 


I . '^CR-Th AN H 

:OMPILE« OPTIONS - NANE^ P.A IN ,OPT = ’’a ,U KECNT = 56 , SIZE =OrCC K , 


SnuKCE, EOCDIC ,NOLISl,NOOI CK ,L C AO , H AP , NOE 01 T , I 0 .NllXK E F 

10e«OUTINE^CCNHAT(NT API,NTAP2,KKcC,NtL,H,N,lX,lV,IH> 

DIHENSiOK 'lX(ijEL't2>,iY(NEL,2l,IhtA,a,M,N) • 

TH'rS^PROCRrH 'KlNdrCONT INGENC V'“PATkICFS IN6 ”ICATTng”AGRcEHENTS 

BETviEFH TWO HAPS IN TERMS CF CLASS LAdfLS AND THE BOUNDARY TYPES. 


IH(<‘.<',Ii J) REFERS TO LOCATIONS WITH CLASS I IN MAP I AND CLASS J 
IN M AP 2. LEFT HALF OF IH GIVES A COUNT OF AGREEMENTS AMD THE 
RIGHT HALF, DISACREEHEMTS.. 

ROMS O F COR RESPOND TO MAJLl^ A ND C LO UMNS T O HAP2. 

h(Tw NUMBERS I, 2 , V, A INDICATE NO BOUNDARY, CHTnGE iN VERTlcit 


- -PATE 7A..259/T,1,,A5,5JI_. 




c 

DIRECTION, CHANGE IN 

HOR IZCMTAL 

DIRECTiJNi 



c 

c 

. TIONS, RESPECTIVELY, 
TYPES OF TRANSITJCNS 

IN HAP 1. 
-Itf„flAEij 

S IMILARLY 

ISN 

poo^ 

c 

11 = 1 



. \'n 

ISN 

bjo's' 

^DOo 


12*2 

NCL^t^NEL^A 


' 



ISA 00 0 7_ 
"iSA 'do'.)3 
.iSN_niG.i. 
ISA jOiJ 
ISA )01I 


ISN 0H2 


IS N JH3 
ISA OTIA 


I-SlL-iil?-. 


ISA JJ17 

.ISL.QJ18. 
i’SN 0Di9 
ISN 012-7 


ISA .0721 


INITIALIZE THE ARRAYS IX AND lY. THE «?REVIGUS« RDM TG RUM 1 IS 


! CQNSIOEkEO IDENTICAL TC ROM I, 

_L£ALL SAKN(NTAP1,IX,NFLAI : , 

' CALL SARN(NTAP2,IY,NELAl 

j_CALL VHQVnx,NEL,IX(l,2)) 

, CALL VHOVI !Y,NEL,iv( 1,2) ) 

• CALL SVSC I ( 1 F.32« M» N .CI 

I 

^lOOP ON RE CORES. _ 

DC 10 I=1,NREC 

"Tqlip'on pi'xels;; ■ "■ 

DC an J=1,N E.L 

! JP=KAXOIl,J~-l) 

T’NFNP’osi’TiVE VALUES OF HAP LABELS ARE' NOT UF’'iTtERE$T 
J_JJFJ-L!LU,.L?i.,Le»0.,DR.lY(J,ia,ULE.C)G.q. TO 2C 

r FINC ROW AND COLUMN NUMBERS IN IH IC BE INCREMENTED. 


1— THE..NATURI LF_TFE BCLNQA^IES IN BOTH _l>lg. HAPS.. 

r K = IXlj,I2) ~ ~ 7" 

' L = IVU.I2) . . 

; Il = l + INTLOGnx(J,ia) .NE. IXIJ, I lll+a’*INTLOGIlx’( j,i 2 ).N£.IX(JP,I 2 i) 
, JJ = ltINTLC) 0 (IYt J.tai.NE. IY(J, 1 1 )) ♦ 2 «'I NTLUG H Y ( J , 1 2 J . N£ . LY I JP , 1 2 ) ) 


i F IND THE INCREMENTS. . . __ 

t] INC = NUMBER OF AGREEMENTS? INC 3 » NOMbER OP Ols AGREE MENT S . ~ “ 
ii.JLNC = INTLPGnx( J,ItI.EQ.IYI.J,I 2 )) „ 

,1 . + INflbGnX(J,ii).EQ.lY(J,Iin ♦ fNTL0GnXUP,'r2).BC.IY(JP“,i2)i 


"op tS" 



'1 







ISN 

0022 


INC3=^3-INC 


ISK 

\su 

0)23 

002A 


IH( n«JJ»K«U= [h(II, JJiK.U * INC 
IHni,JJ+A,K,L MIHHI, JJ«-^,K,LI*INCa 


ISU 

0025 

20 

C 

CONTINUE 


ISU 

0026 

C 

EXCHANGE n AND 12 
IH= 11 


is\ 

I5N 

o)27 

0023 


11^12 

!2«IU 




C 

C 

RE/)C NEXT RECORDS INTO I X( » « 1 2 1 « i V ( «. I 21 


ISK J129 
I Sis l0 31 


IF(UEa,NREC)GU TD 10 
CALL SARN(NTAP1,IX^ U 12) 


ISK 0J32 
JSK CJ3J 
lU 0 j3a" 
rSK 0035 


CALL SARN(NTAP2,IYU, I2l‘fN6LAl 
CCNTINUC 


RETURN 

END 
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I5N J015 • 00 2C 

011^ 20 HRITE<6»3C<^MlH(K»KK>IrJ) iKK^jlII 

iSn u3i7 io comtinUe“ 


TsToyrr 

ISh C319 




FIND AND 
^OR EACH 
OFTLAyS 
DF^EACH_ 
‘do 30““I = 
DO 30 




tINT MATRICES SHGiiING CCLNT5 OF EACH TYPE OF TRANSITION 
UR OF CLASSES.^ THESE^SHOWr FO R ALL JO INT CCCURENCES 
(Tf'yriN HAPS^i/27 fHE MjMyERs’’OF JOINT OCCURENCES'^ 

>E OF TRANSITION IN THE T^G MAPS. 



ISN J-J33 
ISN 

Un jj35 


DO 5C 

Jd ( 6 1 , 3 1C lU H LI* KK I ItJ LrK K =J-i f 

CONTINUE 


ISN O'^3o 
ISN 3337* 
, ISN OJ3 9 
ISN u0 39 


FIND IHl, THE MATRIX OF COUNTS CF TRANSITION TYPES! W ! THOUT REGARD 
TO CLASS LA6ELS> AND IH2» THE MATRIX OF JOINT OCCURENCES OF ClASS 
LABELSIWITH NO REGARD TO TRANSITION TYPES). 

CALL_SVSCn IH1| UjC) 

CALL SVSC n iH2,M^N,C) 

^0^60^ fM 

OG 60 J=1,N 





liK 0 j 68 CAU^ SVSCn INVliH.C) 

is^ 1)36? ■ CALL'TvSC lTlNV2',S,Ci 

ISH o3 7J 00 8 5 1 = 1 , H 

isk'TjTyi oo as j=i7n 

ISh 0)72 [NVU n»lNVl( n*Ih2(I ,J» 



ISA 0079 310 FDKMAT('''*15IB) 

ISA COdO AOO FORHATI/' MATRICES SHOHlNfi COUATS CF ACREEHEMTS AND DISAGREEMENTS 

"■ .FOR eACH"TYPE“0rfRA'NSITI0A»/'“>A~P'siZt = 'I5,* 

ISA 0031 AlO FORMAT!/' MATRICES SHOhING CQLNTS CF EACH TYPE CF TRANSITION'/ 

. /• MAP‘5IZE = 'I5V''8Y'I5f 

ISA 0)82 A20 FORMAT!/' MATRIX SHOWING TOTALS OF EACH TYPE OF TRANSITION'/ 

' MAP SIZE = *I5,' B Y • fS , ' , • I 

ISA_00d3 A3!) FORMAT!//' CUNTIGENCY TABLE') ^ 

TSA OTdA" 500 F0RMAT(7/'“'CLAS1 AUMBck" IN PAP i='i2,' CLASS NUMBER IN MAP 2»'l2) 

ISA C,3d5 600 FORMAT!' T RANSITICN S I P I L A R 1 T I ES= ' 1 7, ' T T0TA L=*I7> *{ PERCENTAGE^' 

F7.2)"' 




I S N C ) ii o 


_ 003 

701 FORHATi* NUHBFR CF POINT BV PCIUT S IM ! L AR I T I E S* M 7 1 • T0TAL«M7f 
« • PFRC£NTA'GtTrp7r2) 

ISK 0337_ B 0^__F D R^fU ^ ( / •_ UVnUCRY UP M AP I : * / U X I 5 I fl ) > _ ^ 

lSK“00'5d 810 PORMAtf/* INVENTORY OF MAP2J«>(lXl5I8n “ ' " 

ISN 0T89 END 



Ci MD Oi^ 


LbVCL ( JUKI 7<t 1 CS/,3tr_ fLUKAI^. H t)ATE_ 76.259/11.46.22 

CCHPRER OPTIQtgS - NAME* HA tN,OPT*02,LlNECNT»56, SIZE-'^'CCK, 

rnllk c E . E e c 0 ic7S 5 CTriTNC ol' c'k Vl c a ti , V! a p , mu t o i t , i o , n o' x r e p — 

2_ J FUNCTICJN lNTi.UC|U _ _ ' 

J LGCiCAL L 

4 1NTLDG=0 

5" • IF(LI 1NTLGC = 1 

I SAJLLOT RETURN 1 ; ! : 
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GTMV/S LCH1120573 HCDILE BAY DATA! 


MATRICEi SHCUINC COUNTS OF AGREFHENTS AND DISAGREEMENTS FCft EACH TYPE i 
HAFTiZts rj2A BY 8b6. ' _ — _ 


C LASS NUMBER IN HAP 1= I CLASS NU MBER IN HAP 2= 1 


7U01 12518 

13950 

4568 

1 

0 

6259 

6979 

9136 

792 

A79 

196 

282 

1 

396 

97 

392 

3C6 

798 

182 

526 

267 

1 

399 

' 36A 

IBA 

309 


77 

92 

14C 

1 

IM 

r>i 

icr 

121 

103 










CLASS NUMBER 

IN MAP 

1 

CLASS NUMBER 

U MAP 

2= 2 





30'B I TIW 

"4 ’8 T9" 

'I "" 

“ 39B5'8 

1291 5" 

■ ■ r3902~' 

■8'56’7“" 

mmm 

96 

TT7 

pfT5 — 

1 

1415 

■ 612 


— jn — 

132 

• ISC 

ITff' 

“173* 

1 - 

17lfc 

ai 

■' • 786 

‘ "’4T3 — 


3^ 

■ 41 

■5 5““ 

r ~ 

■■■■ 435 

" “m" 

— lar" 


CLAiS KUHaeK 

IN MAP 

1 

CLASS NUMBER 

IN MAP 

2*^_3 




2159 

2327 

31A7 

1 

51C39 

^ 12^7^ 

12A33 

6324 

627 

38 

27 5 

73 

1 

1953 

596 

AA5 

_34 7 

775 

316 

31 

5fi 

1 

_ _ 22 73_ 

53C 

692 

353 

527 

5A 

55 

9 

1 

Sq2_ 

_iS3_ 

167 _ 

156 










CLASS hUHtiEK 

IN MAP 

1- 1 

CLASS NUMBER 

IN MAP 

2* 4 



0 

A1 

rsT 

173 

1 

~ 50 40 

*69T" 

lU 3 

■■~6T6 

121 

c 

54 

4 

r 

2o3 

21 

'78 

14 

A6 

23 

6 

1 

1 

lli 

52 

15 

26 


2C 


9 









GTM y/5 LCMn 20572 MOBILE BAY DATA) 

matrices showing counts of each TYPE CF TRANSITION 

• 

- *• ' * 

, . * 

- 

. 


^ ' 


» 1 
1 ‘ 
»t 

_MAF SIZE=^L624 JY 966. 


- 






- 



... .. 

trX 

t*l 

Cl 

CL45S KUI'llfift IN HAP 1= 1 CLASS NUMBEK IN MAP 2= 1 









n 

■'> 



CLASS NUMdER IN MAP 1= 1 CLASSJ.UMBfR IN MAP 2^ 2 



CLASS NUMBER IN MAP 1= I CLASS NUMBER IN MAP 2= 

^ — 1 7.1 r 3 2 T 9 ri! 2TV7I3 


rn6 2 62 




CLASS NUMBER IN HAP 1= I CLASS NUMBER IN MAP 2- A 


1 bJ2 245 A2 3 

44 



128 


7 


6 



C TH V/ 5 LCM(I2C573 HC6ILE BAV DATA) 

HATRU iHCrtlNG TJTALS OF EACH TYPE OF IPANSITIOH 


PAP SI^E= lo24 EV 866. 

■ •■■57r501 i:5S28'3 TlTsTS 77'l71 


T5‘9y5 577F 



NUPriEti CP PlIM OY POINT S IH IL AR IT I ES = <5 7 7 3 39 T01AL = 13A1583 PEKCENTACE® '72.85 


M 

^INVENT CKY__i;F MARJ 

■■l2'lC35 iT379C 522378 A071C2 9677A 10509 


INVENTCBY OF HAF2: 

TA0391 231752 A97515 A0S3A2 91C77 


5511 





9 . COMPABISON OF SUPERVISED CLASSIFICATION MAPS 

9.1 NAME 
COMPSUMP 

9.2 PURPOSE 

To compare two supervised classification maps (or a Ground Truth. Map and 
Supervised Classification Map) and print their joint histogram and the numbers and 
percentages of various types of differences. 

9.3 ■ CALLj[G_SjEQUENCE 


This is a main program with card inputs as follows: 

READ 100, TITLE 
READ 300, NREC, NEL, M, N 
100 FORMAT (72A1) 

300 FORMAT (416) 

where 


TITLE is a title of up to 72 characters (to be printed); 

NREC = Number of records in the two maps; 

NEL = Number of pixels per record; 

M, N are the numbers of classes in maps 1 and 2 . 

9,4 INPUT-OUTPUT 

9.4.1 Input 

The input maps 1 and 2 to this program should be on units 8 and 10 
respectively. They should have NREC records, NEL pixels per record and 4 
bytes per pixel in unformatted FORTRAN readable form. 

9.4.2 Output 

Besides the printout, this program writes difference map on unit 12, with 
NREC unformatted FORTRAN records of NEL pixels each having one byte per 
pixel . 

9.4.3 File Storage 
None 
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9.5 EXITS 


Not applicable 

9.6 USAGE 

This program is in FORTRAN IV and implemented on the IBM 360 with the 
H compiler. The program and the associated subroutines are available as load 
modules on the users’ library. 

9 . 7 EXTERNAL INTERFACES 

The subroutine linkage is indicated below: 


Calling Program 

Programs Called 

COMPSUMP 

SVSCLl 

DCLARE 

DCLARE 

CLAM 

PRTMXP 

SVSCI 

SARN 

VMOV 

CLAM 

SVSCI 

PRTMXP 

MXINX 


9 . 8 PERFORMANCE SPECIFICATIONS 

9.8.1 Storage 

This main program,- presently limited to NEL s 2000 and M*N ^ 200, is 
34932 bytes long, and requires -72K bytes of storage with the external references 
and buffers . 

9.8.2 Execution Time 

The execution time is approximately proportional to the number of pixels 
in the input maps . With NREC =NEL =2000 the program v/ill take approximately 
7 minutes of CPU time on IBM 360/65. 


178 




9.8,3 Kestrictions 


NEL£ 2000, M*N^200 

9.9 METHOD 

This program first sets an M by N matrix LTABL by 

LTABL(I,J) = 'H' forly^J 
LTABL(I, J) = blank for 1= J. 

Tben it calls the routine .DCL^E. This subroutine finds (call to CLAM) and prints 
(call to PRTMXP) the joint histogram between the two maps. 

The next part of DCLARE is used to separate the types of differences between 
the two maps and indicate them by different symbols . The EBCDIC ch^acters 
blank, ’+', and 'H' are used to indicate no difference between the maps, ex- 
terior points, boundary points where the maps are different and interior p oints 

where the maps are different, respectively. The "exterior points" are defined , 

as those where the "class labels" in either of the maps are less than or equal to 
zero. The "boundary points" are those whose class lables are different from 
that of at least one of their four nearest nei^bors (top, left, bottom and right) 
in map 1 . Points which are neither exterior nor boundary points are called 
"interior points". 

These indicators are generated for each of the points in the maps and an 
NREC by NEL pixel sequential data set (unit 12) . The numbers and percentages 
of occurences of these indicators in the output data set are counted and printed 
(except for the exterior points). The percentages'of occurrences are evaluated 
based on all but the exterior points. 

9.10 COMMENTS 

The data set on unit 12 can be - directly to generate a difference 

map. 

9.11 LISTINGS 

The listings of the program and the important subroutines called by it are 
attached at the end of this section. 

9.12 TESTS 

This program has been used in deriving the difference maps and similarity 
measures between several pairs of classification maps and found to work satis- 
factorily. 
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Lp„VFL 21 .S ( JUN 7A \ 


CS/36C FCRTRAN H 


OATH 

CUMPILER QPJl-ONS - NAMT* MA I N tOPT*0 2 ,L I NECNT -56 ,S ! ZE •^'^COK# 

SnURC[!t6DCOICtNnLIST,NOC£CK,LOADtMAP.NOED!Tr IdVnQXREF ' ' 

ISN 'ioa * _„„QJMFNS!QN I X ( ) , l Y ( 20^^: ) . , _ , . 

ISN jOni LOGICAL^*! L Y ( 2:^ }C) » IT ABLI 2C0 1 . , - 

ISN LOGICAL BLANK/^ V»ErCH/»HV 1 

ISK LOCICAK’l TITLE(72) 

.ISt^ Ojjo f^EAp ir'n^TiTL*: : 

ISN '■1)7 PRINT 2:^''. TITLE 

ISh 0U3 READ..?:: tnivEC#NELih,N , . 

ISN L)17 print 4Cr ,NREC ,NEl tM,N 

ISN IM;) _ ^_CAU SVSCLULTA6L|M«NrETCHJ - ! 

ISN sJ'U MN = MINr>(H,N) 

ISN r.112 f)Q 10 I = 1|MN 

ISN ^713 io LTABLni-1 ) )»BLANK 

ISN OU , CALL OCLARECNREC,NEL,H,N,LTA0Lf IXf IY#LV> 

ISN Ojlb STOP 

ISN ons too iORMATaaAU • 

ISN ^U7 200 FOKMATP l'25(/) 3rX72Al) 

ISN n.U 3C0 FORMAT! A )[ 6) ^ ^ : 

ISN >019 400 FDRMATlisX* IMAGE siZE»M5t* 0YM5t*. NUMBERS OF CLASSES IN MAPS I 

. _,..AND.2 ARFM3,* ANDM3t"#M . , 

ISN fj.)27 ‘ END 


00 

o 





,LeVt;L 2K3. (. JUN 7:^) 


05/36C FORTRAN H 


DATE 76,^5?/J6. A8,,57^ 


— 

T C M 

SOURCE, EbCOiC ,NOLIST,NaDkCK,LOAO|HAP,NOEOIT,lD,NgxREF 
SUBROUTINE OCLA R E (NKE C • NEL ,M, N ,LT ABL, I X , ! Y <LY I .. 


' 


t' 

c 

- - -C 

c 

c 

* 

THE PURPOSE OF THIS ROUTINE IS TO FIND AND PRINT: THE JOINT HISTOGRAM , 

BETWEEN TWO MAPS, FACH WITH NKEC LINES AND NEL PIXELS PER LINE. 

THE FIRST HAP SHOULD HAVE R CLASSES OR LESS, AND THE SECUNO, N OR 









c 

c 

LESS. THE ARRAY LTA3L IS USED TO’ SPECIFY ERRONEOUS COMBINATIONS 

OF n A'!'! NllMBPRS. LTABL IS AN INPUT ARRAY TREATED AS AN H BY N MATRIX .. 



' 


c 

_{; 

WITH BLANKS AND H'S. LTAbUIjJI = BLANK INDICATES THAT CLASS 
NUHBERS_I. ANO.J IN HAP 1 AND 2 CORRESPOND. THE INPUT MAPS 1 AND 2 





*^C 

c 

SHOULD BE ON UNITS 8 AND 10 HAVING INTEGER*A NUMBERS C THRU M AND 
0 THRU N RESPECTIVELY (UNFORMATTED FORTRAN). 1_ 





c 

THE OUTPUTS OF THIS ROUTINE ARE: 





. . _c 

n PRINT OF THE JOINT HISTOGRAM. THE NUMBER AMO PERCENTAGE . 




u 

c 

OF CORRECT CL A SS ! F IC AT I DNS , ERRORS AT BOUNDARY POINTS AND 


' 



c 

ERRORS AT INTERIOR POINTS! . .... 




II 

' c 

2) DIFFERENCE HAP ON UNIT 12 CONTAINING COOES IBLANK, .« «■ 




n 

jC 

AND HI WRITTEN AS LOGICAL*! RECORDS. .. 






00 

-HI 


ISN 

ISN 
ISN 
ISN 
rsN 
ISN 
ISN 
ISN 
ISN 
ISN 
I:)N 
ISN 
ISN 
ISN 
ISN 
_„.15N 

IS IV 

ISN 

ISN 

ISN 


ISN 
. ISN 
ISN 
-ISN 
ISN 
-„..-.L5N 
ISN 

ISN 

ISN 

ISN 

ISN 


D ^03 

OJjS 
0 J06 
0007. 
.1 } N 9 
00 J 9 
G)IJ 
0011 
ni2 
GU3 
j014 
^)15 
w T 1 & 
JD17 
•■' J 1 9 
3J21. 
Q721 
Oj23 
C72'* 
0J2«. 


PlHHNStON lX(NEL,3),IVtNEL»,IH(2S6l 

LOGICAL «■! LYlNELlfLTA'JL(H,h»,BOY 

LOGICAL «■! ETCH/*H'/fPLUS/'+'/rBLAAK/' •/ ,00T/*,..'/ 

DATA U/«/ 

CALL_CLAfJ (9,10,NREC,NELtIH,lXfJ_Yi.HtNI. 

REWIND 9 

REWIND 10 . . J. 

PklNT 50: 

500 F0RMATI//5'>X»J0INT HI STOORAH'/ZJ 

CALL PRTHXPlIHfH.NI 

CALL SVSCHIH, 256 ,0 , . 

CALL SARN(9,IX,NEL<-AI 

_.CALL VMQV(lX,NEL,IXll,2n 

00 IP I=l,NRSC 


.IF! I.LT.NREOCALL 
REAOI 10 HY 
_Q0_20. J»1.NEL. 


SARN(8,IXI1,3I,NEL«AI 


IF(!X(J,2J.LE.3.0R,IY(J).LE.O»C0 TO 30 

LYIJl = LTABHIX(J,2»,IYljn 

IFILYIJKEO.DLANKJGO TO 35 

BOY.sIXlJ, 2) .NE. IXIJ.I I.OR. 1XIJ,2I .NE. IX(J,3). 

. .OR. J.GT.I.AND.IXI J>i,2l.NE. IX(J,2l 

. .OR. J,LT.NEL..ANQ..LX( JAl,2).,.NEjJXt J,2». — 


■M27 IF(BOY)LYtJ»-PLUS 

CJ2 J . . GD .10 35 _ 

3)30 • LYU»=OOT 

2331 . .35 IrirLYIJl . 

3332 iKnw)=iHnwui 

J1-33J 2i3 CONTINUE 

v33A CALL VH0V(IXCl,2> .NEL.IXI l,in 

0335 ^CALL. VMQvnXLl,3J ,NEL,IXI1,2» I 

''336 10 HRITE(12»LY 

3J37 N.T.DT.= 0 . . . _ 

0338 DO 65 1=1,256 


10 


. 182 ; 


ISN OJjf ^LV_(ll = l _ . . _ . 

I.Sh. ifCLYtl r. ruiDIJTJCO TO 65 

15N 0162 Nn)T = NTnT + lH(n 

ISN 65 CONTINUE 

. liN 0166 ... FAC=1C0./NT0T . . _ 

C 

C PaifJl. HISTOGRAM ANO..PcRCENTAGE_eCCUPeNClES. . 

C 

ISN 1065. D0..6G I=>.l,25fc _ 

ISN 6166 IFMHm.Eg.OIGO TO 6C 

. . . ISN 1)6d PEBCNT=lHin<‘FAC . ' 

ISN 1569 LYm=I 

ISN.115.1 IF(LY(n.EQ.DOTIGQ TO 60 _ 

ISN 1)52 PRINT ICO, LY(lI,lHm,PERCNT 

. ISN 1)53 ...60 CONTINUE _ _ __ 

ISN 0356 100 F0RMAT(/6rXAl,llXl8,l1XF7.2) 

ISN 0)55 PRINT 200 

ISN 0‘>56 201 FORMAT(//ICX"'.*'= EXTERIOR, « « » NC ERKOR, '♦ 

. J.Y POINT. "H” a ERROR _AT NONBOUNOARY_POINT • I 

ISN 0)57 RETURN . 

ISN. 1153 E.NO_. . . _ . 


BEERODUCBBIL] 

ORIGINAL PAG: 



LevIL 21.a ( JUN.7.A.I CS/36C FORTRAN H j ' DATE. „ 76,25?/.16, $?,(}&. 

• CO MPILER OPTION S - l<AWF» ~ H a IN .OPT gf 2 ,L 1 NEC MT = 56 , i I ZF sOOnc K « „ 

• SOURCC,CPCniC,HOHST,NCnfcCK,LOAD,HAP,NOEOIT,rD,NQXREF 

.. liA OOJ-i SUBROUTINE PRTM XP (I A , H ,M I . . 

"• ISK OOOa DIMENSION I4(M,N> ‘ . 

ISN CALL MXINXI IA,H,N,H1.NJ> . . . - 

I5N 0)DS DO 5D I»1,MI . ' j 

ISH.0J36 SO PRINT ■330»tUtl.«.JI.j£Ji,NJ> - ! 

ISS 0007 3C0 FORMAT! 10X15181 


>■ ISN 0033 RE.TURN 

I5N 0309 END 



00 

. CO- 

I - 


I ^ 



I . 




10 . TEMPOEAL CHANGE DETECTION AND DISPLAY 

10.1 NAME 

CHAEADE 

10.2 PURPOSE 

To provide a difference image coded so as to indicate uniquely the class 
numbers in each of the two input images, and whether each pixel is an interior or 
boundary pixel; and to identify, display, and tabulate particular types of change. 

10.3 CALLING SEQUENCE 

Each of the two tasks described is controlled by a main program, 

10.4 INPUT-OUTPUT 

The inputs are certain parameters and the two data sets containing the 
classification maps . The parameters read are listed for the two tasks . 

.Difference Coding; 

NSL: Number of scan lines in the map, 

IS; number of samples per line, 

M, N: maximum numbers on the two input map data sets, 

ITAB: ’:look-up table to reassign the class numbers in one of the maps, 
IFLG: an integer which was added to the class numbers at all 
boundary points (see program FLGBDRYS) , 

The first two of these parameters are read from cards; the remainder have been 
written on external storage by the program MXSMLRTY. 

Change Detection: 

NREC: number of records in the data, 

NEL: number of elements per record, 

M; number of classes in the maps, 

IBUF: an array containing codes specifying the types of changes 
to be detected and displayed. 

The output is the difference coded map, a printer plot and external storage 
file depicting the changes, and an inventory of the changed pixels. 
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10.5 


EXITS 


Not Applicable 

10.6 USAGE 

The programs are written in FORTRAN and were compiled using the E 
compiler on the IBM 360. They are on the user library in executable form. 

10 . 7 EXTERNAL INTERFACES 

Certain other subroutines are called for generating the change depiction 
image on the printer, for I/O, and for vector manipulation. These are also 
available on the user library, and are the following subroutines: 

PLTPIX 

SARN 

SAWN 

SVSCI 

VMOV 

10 . 8 PERFORMANCE SPECIFICATIONS 
Storage: 

Difference coding - 42K byte s (K=t1024) _ 

Change detection - 56K b 3 rtes 

Execution Time: 

For map si25es 1624 records by 866 samples (ij 406,~38Tpixei^r. 

Difference coding - 1 minute 
Change detection - 1 minute 

I/O Load: 

The difference coded image is placed on external storage for passing 
between steps. The change detection image is written on external storage 
for use in plotting on the printer . 

10.9 METHOD 

The two classification maps are read, class numbers are reassigned by 
table look-up if desired, and the difference coded map is written on external 
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storage. The user's request array is read to determine the type of change to be 
considered. The difference coded map is read in order to produce an inventory 
of the changes as well as a change depiction image. This output image is written 
on external storage and displayed as a printer plot. 

10.10 COMMENTS 


This section is devoted to a guide in the coding of the user's request 

array. 

The user is expected to input his request of change detection and display 
through a data card (using the format (2613)) which is read into the array IBUF 
by the system. Accordingly, the first entry into the array, IBUP (1), should be 
set to indicate the user's region of interest using the following code. 

IBUF (1) = 1: all re^ons (pixels) in the scene 

= 2: interior regions (pixels) only 

= 3; boundary pixels only 


The second entry, -IBUF (2),is set to denote the user's desire as to the type 
or category of change to be depicted using the ensuing code. 


IBUF (2) = 1: 

= 2 : 
= 3: 

= 4: 

= 5: 


each change from each class (in the set to be defined 
for Map 1) to a corresponding individual class (in the 
set to be defined for Map 2) . 

change from each class (in the set to be defined for 
Map 1) to every class (other than itself in the set to 
be defined for Map 2) . 

change from the union of classes (in the set to be de- 
fined for Map 1) to the union of classes (in the set to 
be defined for Map 2) excepting of course those pixels 
that remain unchanged. 

change from the union of classes (in the set to be de- 
fined for Map 1) to individual classes (in the set to be 
defined for Map 2) disregarding of course the pixels 
that were originally in the corresponding class . 

change from individual classes (in the set to be defined 
for Map 1) to union of classes (in the set to be defined 
for Map 2) excepting those pixels that remain in the 
same class. 
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NC1=IBUF (3), the third entry defines the number of classes to be con- 
sidered in Map 1 with IBUF (4) through IBUF O^Cl+3) specifying the actual class 
numbers in the set of classes. However, if IBUF (3) is set to zero, then this 
- denotes-that a sequential set of classes is to be considered and in that case 
IBUF (4) defines the number of classes to be defined internally in sequential order 
for Map 1. Similarly, the next entry NC2=IBUF (NEXT) (where 

NEXT = 4 + NC 1, IBUF (3) 0 

= 5, -0 

defines the number of classes to be considered in map 2) with IBUF (NEXT + 1) to 
IBUF (NEXT+NC2) the actual class numbers in the second set. As before, if 
IBUF (NEXT) - 0, then this denotes that a sequential set of classes is to be con- 
sidered in Map 2 and in that case IBUF (NEXT+M) defines the number of classes 
to be defined internally in sequential order for Map 2. 
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Typical Examples to illustrate the Code for input of User's request: 

1. User's request: Depiction of the Union of Changes occurring in the 
entire scene from all 9 classes in Map 1 to all 9 classes in Map 2 . 

The code -will be 


/ 


denotes total 
changes 
(interior and 
boundary) 


denote a set of sequentially 

numbered classes in Map 1 and 2 




denotes changes 
from union of 
classes specified 
for Map 1 to union 
of classes specified 
for Map 2 


number of 
classes of 
Map 1 and 2 



2. User's request; Depiction of changes in the interior of Class 1 (Map 
to Class 2 and 3 (Map 2) separately. 


1 ) 


The code is 


number of classes 
of map -1 - 

\ 

2 . 1 ■ . 


set of classes 
- for map 1 

1 2 


denotes 

denotes changes 

number of 

changes 

from each of 

classes of 

in the 

classes specified 

Map 2 

interior 

for Map 1 to every 


regions 

class specified for 


only 

Map 2 



set of 
classes 
for Map 
2 


REPRODUCIBILmOFm 

miCINAL PAGE IS POOR 
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3 . User’s request; Depicftion of changes in the boundaries from Class 1 to 
Class 2 and Class 2 to Class 1 separately. 

The code for this case is: 


number of classes 


. . 3 .. . 



denotes 
changes in 
boundartes 
only 



from a class 
specificied for 
Map 1 to a 
corresponding 
class specified 
for Map 2 


set of classes 
for Map 1 

1 - « • 2 • 


«2* «2» «X 


number set of classes 
of . - for Map 2 

classes 
of- Map- 
•2 


4. User’s request: Depiction of changes occuring in the interior of classes 
1 through 5 taken together into classes 3 and 6 individually. 


The code for this case is: 


. . 2 

/- 

denotes - • • 
changes in 
interior 
regions 
only 


denotes a sequential set of classes 

set of classes for Map 1 for Map 2 



denotes changes 
from union of 
classes specified 
for Map 1 to 
individually specified 
classes of Map 2 


number 

of 

classes in- 
this 

sequential 

set 


number of 
classes ' 
of Map 2 
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5. 


User’s request: Depiction of total changes from Class 1 and 3 individually 
to Classes 3, 4, and 5 taken together. 

The code for this case is: 


/ 

denotes 

total 

changes 


number of 
classes for 
map 1 

. 5 . . 2 . . 


set of classes 
for Map 1 

i • 

J. * ♦ o . 



2 . . 4 


5 


denotes changes 
from individually 
specified classes 
of Map 1 to union 
of classes specified 
for Map 2 


number .• T set of classes for 
of Map 2 

classes 
of Map 
2 
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10.11 LISTING 


- 

1 n C O 1 U A_\ it; U ♦ l ^ j. V f 

DATA IN, JNrKN/8 ,9,10/ 
READfS-.in.) NSL.IS 

♦ 1 1 H U t > — 


10 

F0RHAT(2I6) 

R BADm M . N , (I T A8 ( T ) ,^1 =! , N ) 

,LFLG 


M = M - IFLG 
WRTTE(KN) M 

IOC 

CALL 0VRMSK(TN,JN,KN, 

FORHATdOX, 'DIFFERENCE MAP 

IX, lY, TT4B,NSL,IS,M) 

CODING HAS BEEN COMPLETED WITH M = ' , 1 3.) 


WRITE(6,1C0) M 
STOP 


END 


SUBROUTINE .0VRMSK{IN,jN,KN, IXtIY, ITA8,N5LtlS,N) 

DIHENSIO.N, I.XnS } ,IY{ IS) , ITA8( 1) 

C TO GENERATE A MASK WHICH INDICATES WETHER A GIVEN PIXEL IS INTERIOR 
C (OR BOUNDARY) AND THE CLASSES TO WHICH IT IS ASSIGNED IN THE TWO .MAPS 
C ON UNITS TN,JN. THE OUTPUT APPEARS ON KN. TTAB IS A TABLE INDICATING 
C the assignment of CLASSES IN MAP 2 TO MAP l^IFLG CONSTANT ADDED TO 

C CLASS NUMBERS IN HAP 1 TO INDICATE THAT A GIVEN PI X E L IS ON THE 

C BOUNDARY. • - ... 

DO II I ^ 1 ,NSL ~ 

PEADIIN) IX 
READ(JN) lY 
DO 12 J = 1,IS 

IFHXI J ) . EQ.O) GO TO 12 ; 

I F{ IY(J l.LE.G) GO TO 13 

_C IX(J).EQ.(? INDICATES THE POINT IS OUTSIDE THE REGION Q.= INTEREST 

IX(J) = UX(J)-1)*M + ITABUYU)) 

GO TO 12 . • 

13 IX(J) = 

12 CONTINUE 

WRITE (KN) IX 

RETURN 

END 


reproducibility of the 

GBI05NAL PAGE IS POOR 
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0 T m'f'm R TOW I X n non r. I H I 2 5 6 ) , I S ET ( 1 0 ) . J S.E T { 10) . K5 E T ( I-' ) . I BU F ( ?6J 


INTEGER TABLE(250) 
LOG TO at =!=i 1 xM non) 


CALL CHANGE! IX » L X ♦ IS ET » J SET ,K SET , IH ,I BUF ,TA8LE ) 

RTOP 


END 



SUBROUT INF VN A Y'.S ( -TX « N ) 

DIMENSION IX(N) 

nn TO T = 1 ,N 

10 ix( n = I 

: R ETI JRN 

END 
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DIMENSION IX{1000),IH(256),ISET(1G),JSET11C) ,KSETU0) ,IBUF(26) 

LOGICAL I.xnnon) 

INTEGER TABLEtZSC) 

REAnfS.in )._NIREC fNFI 

READ(20) M 

-■■9 0 Q FORMATHHl. ! 0 y. »THF TNPIIT CHANGE OEPTCTTON RFOUF.ST COnE TS;»«? 6 T 3 > 

800 F0RHAT(26I3) 

10 FORMATf^Tfel 


READ(5,3Ca,EN0=l) IBUF 

PEFTNF SFTS OF Cl TO Rg CONSTOFi 

IF(IBUF(3).EQ.O) GO TO 30 
NCI = TRlJFfTl 


CALL VMQV {IBUF{ A 1 ,NC 1 , ISET) 


GO TO 40 

3n NCI - TRUFI^l 

CALL VNATS( ISET ,NC 1 .) 


40 

IF{ IBUF( NEXT) .EC.O) GO TO 5C 
NC2 ^ IBUF(NEXT) 

CALL VHOVnBUF(NEXT+l) ,NC2,J5ET) 
GO TO 60 

50 

NC2 = I6UF(NEXT+1) 
CALL VNALS-I JSET,NC2) 

60 

00 15 L = 1,26 

1 = ?7-L 


IF(IBUF{ I ).N£.0 ) GD TO 16 

IS 

rnWT TMUF 

16 

WRITE{6,9C0) ( IBUF( J) ,J = 1 ,n 

Kwlira 

TSETH 1, T=1 .NC 1 AND J SFT f T) . 2= 1 . NC ? APE THE CiASS NUMBERS USED 


KCl = NCI 



DO 30 L = ItNREC 

CAfl ARM f ?fl . TY .NFI *4 1 


DO 85 LL = 1 ,NEL 

LY fLL‘1 = TABI F < TYM I 1 + 11 



CALL PRTLBLl IBU F , ISET , J5 ET , KSET , NC 1 ,NC2 , I f J,IHJ 







201 

2n? 


2Q3 
ri 


301 

303 


im 




DIMENS ION I8UE{ 1 J f ISETINCn »J5ET( NC2) 
1 J .K'JFT ( 1 1 - 


FORHfliT ( IQX, ’CHANGE DEPICTION MAP SHOWING THE ENTIRE PIXEL SET OF’) 

,F " ^ 


FQRMATdOX, ’CHANGE DEPICTION MAP SHOWING TH 

FfiRMATf / / . 1 nv • wn . flF PTyFl<;« 


FORMAT! 1H + ,63X» ‘MAP 1 - C L A SS ( E S ) : ’ , 1 5 1 3 ) 

■ FORMATfiry . TI^. FX . »nvPRPft T NT» . 7X. «TQ._CL A S S ( E S ) :« .1513 
FORMAT! lOX, ’( EXCLUDING OFCOURSE THOSE PIXELS THAT REMAIN IN THEIR 
.nRTr,7NAt n &<;«:f ) 


FORMAT! IQ X, I13,<?X,'I’tl2X, ’NONE') 

MAT! iny. T 1 PY . » <« » .11 X . ‘Tn Ail FTHF: 


FORMAT!//, IGX, * THE REST OF THE SCEIvE OF INTEREST IS SHOWN AS BLANK 

F THE SCENE DEPICTED BY PLUS SIGN’) 

IBl = IBUF! 1) 


ifsiiiai 


GO TO (10,20,30,181 



GO TO 99 

T 301 . 


PRINT 2GA 
ET 


GO TO 99 
PRINT 3’--l.(TSF 
PRINT 204 
PRINT 305. TH!’ 


£^_i5.EX!Jn PRINT 302. I H (2 5 6 


SFT! U .L= l.NCl) 

) 


I ?nfc ) . { JSFTf L J .1 = 1 . 


PRINT 303 
GO TH 99 


LK = G 

00 93 L~ l.NCl 


IFdSETIU.EG.JSETU) ) GO TO 93 
LK = LKj- 1 


KSETCLK) = ISET(L) 
. CONTINUE 


PRINT 301,!K5ET(L1 ,L=1,LK)' 

PRTMT. 204 

PRINT 3G2,IH! 2561 ,JS£T( J) 

GO TO OQ 


LK = n 

PR INT 301 . ISETT I ) 


PRINT 204 
DO 941= 1.NC2 


r( JSET!L ) .EG. ISETd ) ) GO TO 94 



ORIGINAL PAGE IS POOR 
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t K = I K » 1 

KSETILK) =J5ET(L) 

94 CONTINU E 

P R*I NT-30-2TTHr2 5'6 ) rCKTETC L ) , L = 1 , LK ) 

2 PRINT 3n4.TH(31 ) 

IF{IH(46).NE.O) PRINT 305,I.H(46) 

PRINT 

ITOT = + IH{256I+IH{ SI > 

PTT ^ 7H{2q6)»10ri.n/TTQT 

WRITE(6,7C0) PCT 

mu F0RM_AT(10X, ’CHANGE OF INTEREST: * .F_5.2 .7 ^« ) 

•RETURN 

£Kj] : 
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SUBRnUTINE CREATF( I . J , I BU F , TA 3L E . T SET ♦ J S E T . NCI . NC 2 .M ) 

dimension I8UF( 2) ,TASLE(1 > ,ISET(NC1), JSET(NC2) 

INT = f,FR TiBI F 

C ■■ DIMENSION TiBLElMC33=MCl’!=MC2*2 +1) WHERE MCI AMD MC2 ARE' THE 

_2 NUHBJ=R QF n TN map 1 AND ? 

MS = {M+-U<'M • 

MM = MS + M»M 

!31 = T8UF{ n . 

132 = TRUF(2r 

TABLE (1) = 6C 

CALL. SVSrn TARI EIPWMM.r) 

12 = I 

I3_= I 

J2 = J . 

J3 = J 

GO TO (2C ,3'',^'} ,A0»60) ,132 ' 

2H .17 - T 

J3 = I 

GO TO 3D 

4P 12 = 1 

T3 = NCI 

I82.EQ.A) GO TO 30 

^0 .12 = 1 

J3 = NC2 

30 CONTI NUF 

DO 12 II ^ 12,13 

L = ( rSET ( TD-l }»M +1 

00 11 J1 = 1,M 

K = I + n 

I F{I5ET{ I ! ) .EQ. J1 ) GG TO 15 

IF(IB1.LT.3) TABLF(K) = 45 

IF{IB1.NE.2) TA3LF(K + MS) == A5 

G_C_T0_1_L ^ 

15 CONTINUE 

I P f 181 .LT . 3 ) TABLFf K ) ~ 3^ 

IF(IS1.NE.2) TA0LE(K + M5) = 8C- 

n CONTINUE 

00 I'O J1 - J2,J3 

K = L + J SET {.ill . 

TFIISETd n .pg. JSET( Jl) ) GG TO 10 

IFHB1.lt. 3> TABLF(K) = 255 

IF(I81.NE.2) TA3LE(K+M5) = 255 

CONTINUE 

12 CONTINUE 

RETURN 

END 
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